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Abstract 

This  study  was  motivated  by  a  need  to  develop  a  reliable  method  of  predicting 
the  agility  characteristics  of  various  aircraft.  To  fully  investigate  the  agility  of  an 
aircraft,  maneuvers  which  push  the  limits  of  an  aircraft’s  maneuvering  capabilities 
must  be  simulated.  In  these  cases,  classic  trajectory  optimization  techniques  either 
require  too  many  assumptions  for  a  realistic  solution  or  require  a  good  guess  of  the 
final  solution  before  the  problem  is  even  attempted.  This  study  investigated  both  the 
utility  of  pseudospectral  optimization  methods  for  robust  trajectory  optimization  as 
well  as  the  potential  for  demonstrating  differences  in  aircraft  agility  characteristics  of 
several  specific  maneuvers. 

Building  off  of  a  pseudospectral  optimization  software  package  named  DIDO, 
a  robust  maneuver  definition  and  trajectory  optimization  system  was  developed  to 
simulate  various  maneuvers  specifically  designed  to  demonstrate  aircraft  maneuver¬ 
ing  limits.  This  system  was  used  to  optimize  the  trajectories  of  three  variations  of 
a  baseline  F-16  mathematical  model  developed  to  simulate  important  differences  in 
aircraft  agility  characteristics.  Initial  results  showed  significant  instabilities  in  the 
interface  between  the  mathematical  model  and  the  optimization  scheme.  These  in¬ 
stabilities  were  mitigated  through  modifications  of  the  system’s  cost  function  and  the 
resulting  trajectories  demonstrated  the  relative  advantages  which  can  be  created  by 
subtle  differences  in  aircraft  designs. 

Future  work  in  this  area  should  include  further  refinement  of  the  driving  cost 
function  and  creation  of  a  graphical  user  interface  to  simplify  the  maneuver  definition 
process.  The  resulting  system  could  be  highly  useful  in  other  trajectory  optimiza¬ 
tion  research  as  well  as  non-related  areas  such  as  accident  investigation  and  reverse 
engineering. 
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I.  Introduction 

1.1  Aircraft  Performance  Comparisons 

For  reasons  too  numerous  to  fully  address,  pilots  and  aircraft  designers  have  been 
attempting  to  compare  the  characteristics  of  different  aircraft  since  long  before  aircraft 
were  even  viable  modes  of  transportation.  The  characteristics  used  for  comparing  two 
aircraft  vary  greatly  depending  on  the  intended  audience  and  the  specific  reason  for 
the  comparison  study.  These  characteristics  can  range  from  something  as  simple  as  the 
weight  of  the  aircraft  to  something  as  abstract  as  an  aircraft’s  combat  effectiveness. 
In  the  world  of  high-performance  aerobatic  and  military  aircraft,  the  comparison  of 
various  aircraft  usually  revolve  around  an  aircraft’s  performance  characteristics. 

Initially,  the  majority  of  these  comparisons  were  of  basic  performance  measures 
of  merit.  These  included  classic  measures  such  as  climb  rates,  turn  rates,  speed, 
acceleration,  and  range  [8] .  Though  most  of  these  parameters  are  based  on  an  assumed 
steady  state,  they  are  well  understood,  are  fairly  easy  to  determine,  and  adequately 
describe  the  flight  regimes  of  most  aircraft. 

These  comparison  tools  break  down  when  attempting  to  compare  the  highly 
transient  motion  of  high  performance  aircraft,  particularly  when  involved  in  close 
range  combat.  The  advent  of  energy-maneuverability  comparisons  partially  alleviated 
the  measure  of  merit  deficiency,  but  studies  in  the  early  1990’s  began  to  suggest 
that  these  still  failed  to  quantify  the  ability  of  an  aircraft  to  rapidly  and  accurately 
transition  from  one  flight  condition  to  another,  often  defined  as  an  aircraft’s  agility 
[10]. 
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1.2  Aircraft  Agility 

Whether  sitting  in  a  cockpit  engaged  in  air-to-air  combat,  taking  aim  at  an 
enemy  aircraft  with  a  shoulder  launched  surface-to-air  missile,  or  attempting  to  coor¬ 
dinate  the  impatient  traffic  at  a  congested  civilian  airport,  how  quickly  and  precisely 
a  pilot  is  capable  of  changing  the  state  of  their  aircraft  is  a  very  important  piece  of 
information  for  everyone  involved.  The  ability  to  rapidly  change  an  aircraft’s  state, 
known  as  an  aircraft’s  agility,  has  been  the  focus  of  numerous  studies  and  research 
efforts  in  the  last  two  decades.  Though  a  great  deal  of  effort  has  been  spent  develop¬ 
ing  a  common  definition  of  an  aircraft’s  agility  and  a  set  of  metrics  to  quantify  that 
agility,  there  still  exists  a  need  for  a  method  of  accurately  predicting  and  comparing 
the  dynamic  performance  and  agility  of  various  aircraft. 

Assuming  that  the  aircraft  to  be  compared  are  available  and  that  money  is  not  a 
concern,  a  series  of  flight  tests  could  be  developed  to  run  various  aircraft  through  the 
same  maneuvers  to  determine  which  aircraft  can  achieve  a  specific  maneuver  faster 
than  the  others.  The  real  problem  at  this  point  is  that  it  is  perfectly  feasible  for  the 
best  way  for  one  aircraft  to  achieve  a  certain  maneuver  to  be  drastically  different  from 
the  best  method  for  another  aircraft.  A  simple  example  of  this  phenomenon  would  be 
to  consider  two  propeller  aircraft  which  are  exactly  the  same  except  for  the  direction 
which  the  propeller  spins.  If  the  target  maneuver  is  a  360°  roll  in  minimum  time, 
the  aircraft  with  a  standard  propeller  configuration  would  roll  clockwise  to  benefit 
from  the  torque  from  the  engine.  A  clockwise  roll  for  the  non-standard  aircraft  would 
actually  be  fighting  the  torque  from  the  engine  and  would  result  in  a  larger  time 
required  to  complete  the  maneuver.  If,  on  the  other  hand,  the  non-standard  aircraft 
were  to  roll  in  a  counter-clockwise  direction  it  would  complete  the  maneuver  in  exactly 
the  same  amount  of  time  as  the  standard  configuration  aircraft. 

Unfortunately,  funding  and  time  are  almost  always  severely  limited  and  depend¬ 
ing  on  the  application,  the  aircraft  in  question  are  most  likely  unavailable  or  even  still 
on  the  drawing  board.  These  reasons,  among  others,  necessitate  the  use  of  computer 
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simulations  in  predicting  aircraft  agility  characteristics.  The  most  common  method 
for  accomplishing  this  is  through  the  use  of  trajectory  optimization. 

1.3  Trajectory  Optimization 

The  basic  goal  of  trajectory  optimization  is  to  determine  the  “best”  way  for  an 
object  to  move  from  Point  A  to  Point  B  while  both  minimizing  some  performance 
index  and  adhering  to  the  object’s  basic  equations  of  motion.  The  complexity  of 
trajectory  optimization  arises  when  one  considers  that  Point  A  and  Point  B  are  not 
necessarily  points  in  space,  but  actually  states  which  are  each  defined  by  a  set  of 
variables  and  that  the  performance  index  is  potentially  a  complex  function  of  those 
same  state  variables  as  well  as  the  system’s  time  and  control  variables.  Furthermore, 
one  could  define  a  set  of  boundaries  which  define  limits  for  any  or  all  of  the  state,  con¬ 
trol,  and  time  variables.  These  boimdaries  may  be  system  limits  such  as  a  minimum 
speed,  environmental  constraints  such  as  an  obstacle,  or  performance  limits  such  as 
a  maximum  amount  of  time  allowed  to  complete  the  maneuver. 

Though  the  trajectory  optimization  problem  is  well  known  and  is  often  used 
in  the  field  of  aircraft  agility  predictions,  classic  trajectory  optimization  techniques 
either  require  significant  assumptions  or  a  good  guess  of  the  final  solution  before  the 
problem  can  even  be  attempted.  The  assumptions  required  for  classic  optimization 
methods  rule  out  the  ability  to  optimize  a  full  aircraft  mathematical  model  and  the 
guess  restrictions  place  severe  limitations  of  the  ability  to  detect  drastic  maneuver 
differences  between  aircraft.  Recent  advancements  in  the  field  of  pseudospectral  opti¬ 
mization  methods  now  allow  for  the  optimization  of  full  6-Degree-of- Freedom  models 
with  minimal  assumption  and  guess  requirements. 

1.4  Research  Objectives 

The  focus  of  this  effort  was  aimed  at  developing  a  method  of  controlling  an 
aircraft  independent  flight  simulator  for  use  as  a  tool  for  comparing  the  flight  char¬ 
acteristics  of  various  aircraft.  Previous  attempts  at  providing  external  controls  had 
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resulted  in  only  simulating  fairly  benign  maneuvers.  This  study  investigated  both  the 
utility  of  pseudospectral  optimization  methods  for  robust  trajectory  optimization  as 
well  as  the  potential  for  demonstrating  differences  in  aircraft  agility  characteristics 
of  several  specific  maneuvers.  The  goal  of  this  research  is  to  develop  a  trajectory 
optimization  system  which  will  allow  a  user  to  investigate  and  compare  the  agility 
characteristics  of  various  aircraft  by  simulating  a  wide  range  of  maneuvers.  Once  opti¬ 
mal  trajectories  have  been  determined,  these  results  can  be  used  as  the  control  inputs 
to  a  flight  simulator  model  for  visualization  and  more  detailed  analysis  purposes. 

1.5  Thesis  Overview 

In  Chapter  2,  further  discussion  on  the  topic  of  aircraft  agility  is  provided  along 
with  a  theoretical  development  of  the  trajectory  optimization  problem  and  the  basic 
aircraft  equations  of  motion.  Chapter  3  details  the  specific  aircraft  mathematical 
model  and  trajectory  optimization  methods  used  in  this  research  as  well  as  an  overview 
of  the  specific  maneuvers  which  were  simulated.  The  results  and  subsequent  analysis 
of  the  optimization  runs  are  presented  in  Chapter  4.  Finally,  Chapter  5  contains  the 
conclusions  and  recommendations  for  future  work  which  resulted  from  this  research. 
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II.  Theoretical  Development 

2.1  Previous  Research  and  Motivation 

When  attempting  to  compare  the  performance  capabilities  of  different  aircraft, 
the  parameters  used  generally  fall  into  two  categories:  point  or  integral  [8].  Point 
parameters  are  those  which  are  valid  at  a  specific  point  in  time  and  do  not  take  into 
account  what  occurs  at  any  other  point  in  time,  whereas  integral  parameters  take 
into  account  changes  in  the  aircraft’s  state  over  time.  Examples  of  point  parameters 
include  most  of  the  classic  aircraft  performance  terms  including  stall  speed,  maximum 
turn  rate,  maximum  speed,  and  minimum  drag  speed.  Of  the  classic  performance 
terms,  the  notable  exceptions  to  the  point  parameter  generalization  are  the  range 
and  endurance  terms.  To  determine  an  aircraft’s  maximum  range  and  endurance,  one 
must  integrate  the  aircraft’s  state  over  time. 

Each  of  these  classic  performance  parameters  have  one  basic  assumption  in 
common:  steady  state  flight  conditions.  Since  flight  profiles  for  the  vast  majority  of 
aircraft  are  generally  dominated  with  large  portions  of  steady  state  flight,  it  clearly 
makes  sense  to  base  the  fundamental  comparison  parameters  on  that  assumption.  On 
the  other  hand,  as  combat  aircraft  technologies  have  advanced  and  aircraft  capabilities 
have  increased,  the  need  to  quantify  an  aircraft’s  performance  characteristics  in  non¬ 
steady  flight  regimes  has  also  increased.  Efforts  to  characterize  an  aircraft’s  “agility” 
resulted  from  this  need  to  quantify  non-steady  flight  characteristics. 

Aircraft  “agility”  is  both  a  widely  used  and  widely  debated  term  due  to  the  sim¬ 
ple  fact  that  researchers  have  failed  to  agree  on  one  universally  accepted  definition. 
Proposed  definitions  for  “agility”  have  ranged  from  one  extreme  of  sticking  to  a  dic¬ 
tionary  definition  of  agility  [10]  to  simply  observing  how  successful  an  aircraft  is  in 
combat  [27]  at  the  other  extreme.  The  one  unifying  factor  among  all  of  the  various  def¬ 
initions  is  that  they  all  encompass  some  conglomeration  of  dynamics,  maneuverability, 
performance,  and  flying  qualities  [10].  For  the  purposes  of  this  research,  agility  is  de¬ 
fined  as  the  ability  of  an  aircraft  to  rapidly  and  accurately  transition  from  one  state 
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to  another.  This  definition  is  an  amalgamation  of  several  proposed  definitions  [8,17], 
and  will  be  used  throughout  the  rest  of  this  paper. 

Research  in  the  area  of  aircraft  agility  has  mostly  focused  on  two  distinct  ar¬ 
eas:  developing  a  set  of  agility  metrics  and  developing  methods  of  determining  those 
metrics.  While  the  various  proposed  agility  metrics  are  not  the  focus  of  this  effort, 
they  should  be  mentioned  in  the  context  that,  just  like  the  various  definitions  for 
agility,  the  proposed  metrics  are  just  as  numerous  and  wide  ranging  [17].  Many  of 
the  proposed  metrics  are  based  on  how  quickly  an  aircraft  can  accomplish  a  specific 
maneuver,  with  maneuvers  ranging  from  90°  bank  angle  captures  to  maximum  accel¬ 
erations  turns  followed  by  regaining  lost  energy.  Other  metrics  are  combinations  of 
energy  states  and  aircraft  physical  properties  such  as  wing  area  and  load  factors,  while 
others  stills  are  based  on  relative  pointing  positions  between  two  adversary  aircraft  in 
a  dogfight.  The  main  point  to  take  away  from  the  wide  variety  of  proposed  metrics  is 
that,  until  all  interested  parties  decide  on  a  standardized  set  of  metrics,  any  method 
of  determining  those  metrics  must  be  able  to  cover  the  gamut  of  metrics  or  potentially 
become  useless  if  its  specific  metrics  are  not  chosen. 

As  is  the  case  when  attempting  to  determine  most  aircraft  characteristics,  there 
are  basically  two  proposed  methods  for  determining  agility  metrics:  modeling  &  sim¬ 
ulation;  and  flight  testing.  While  a  great  deal  of  effort  has  gone  into  developing 
techniques  for  determining  agility  metrics  through  flight  test,  the  cost  and  complex¬ 
ity  of  obtaining  repeatable  flight  test  data  for  one  aircraft,  let  alone  numerous  aircraft, 
prohibits  the  estimation  of  aircraft  agility,  especially  if  an  aircraft  in  question  is  not 
readily  available  or  not  yet  even  fully  designed.  As  is  usually  the  case,  a  simple 
cost/benefit  analysis  points  towards  modeling  &  simulation  as  the  best  option  for 
rapidly  determining  the  agility  characteristics  of  a  wide  variety  of  aircraft. 

As  many  of  the  proposed  agility  metrics  deal  with  how  quickly  an  aircraft  can 
perform  a  specified  maneuver,  many  of  the  simulation  methods  focus  on  determining 
optimal  maneuvers  and  limits  to  an  aircraft’s  maneuvering  capabilities.  To  this  end, 


6 


trajectory  optimization  techniques  have  been  used  by  numerous  researchers  working 
towards  methods  for  determining  aircraft  agility.  In  one  such  endeavor,  Bocvarov 
[5]  investigated  time-optimal  reorientation  maneuvers  and  the  benefit  which  thrust 
vectoring  control  could  provide  to  those  maneuvers.  The  main  concern  with  this 
research,  as  with  much  of  the  previous  research,  is  that  the  assumptions  which  were 
made  in  an  effort  to  make  the  problem  more  feasible,  invalidate  that  method  for  a  wide 
variety  of  situations.  In  this  instance,  among  other  simplifying  assumptions,  Bocvarov 
neglected  translational  motion  of  the  aircraft  and  only  looked  at  rotational  motion. 
This  assumption  was  based  on  previous  research  which  suggested  that,  during  rapid 
maneuvers,  an  aircraft’s  center  of  gravity  was  relatively  stationary  in  comparison  with 
the  changes  that  the  aircraft’s  attitude  underwent. 

In  another  effort,  the  authors  address  the  issue  that  although  assumptions  of  a 
point  mass  model  facilitate  getting  a  solution  for  the  trajectory  optimization  problem, 
those  same  assumptions  invalidate  the  solution  results  [9] .  Since  the  attitude  dynamics 
of  a  point-mass  model  do  not  adhere  to  Newtonian  Mechanics,  it  is  a  poor  choice 
for  a  trajectory  optimization  routine  because,  unlike  a  real  aircraft,  the  model  is 
capable  of  instantaneously  changing  it’s  attitude,  which  makes  a  study  of  agility 
a  trivial  endeavor.  Instead,  the  authors  modify  a  point-mass  model  to  take  into 
account  the  fact  that  the  forces  acting  on  the  vehicle  and  thus  the  attitude  cannot 
be  changed  instantaneously.  The  results  are  promising,  but  again  the  results  are  not 
very  representative  of  an  actual  aircraft,  and  if  necessary  could  not  be  used  to  actually 
control  a  non-linear  simulation. 

In  a  followup  effort  to  his  earlier  work,  Bocvarov  [4]  presents  the  results  of 
optimizing  two  heading  reversal  maneuvers  with  a  model  based  on  the  F/A-18  High 
Angle-of-Attack  Research  Vehicle.  The  results  from  these  two  maneuvers,  which 
will  be  addressed  again  later,  depict  the  stark  difference  between  the  results  of  two 
fairly  similar  maneuvers  when  seeking  an  optimum  trajectory.  In  the  first  maneuver, 
where  the  aircraft  is  simply  attempting  to  reverse  it’s  heading,  the  results  show  the 
aircraft  should  roll  inverted  then  pull  through  until  heading  the  opposite  direction. 
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The  second  maneuver  similarly  requires  the  aircraft  to  reverse  it’s  heading,  but  also 
stipulates  that  the  aircraft  must  end  the  maneuver  at  it’s  initial  position  and  velocity. 


The  results  from  this  run  show  the  aircraft  performing  a  climb  and  descent  along  an 
arc  which  returns  to  it’s  original  position.  Once  more,  the  results  from  these  methods 


are  quite  revealing,  but  this  group  also  assumed  a  near  point-mass  model  and  the 
control  history  results  are  nearly  meaningless  if  one  wants  to  control  a  full  6-Degree 
of  Freedom(DOF)  model  through  the  same  maneuvers. 

Each  of  these  efforts  have  shown  that,  given  a  certain  number  of  simplifying 
assumptions,  one  can  gain  a  decent  understanding  of  an  aircraft’s  inherent  agility 
through  trajectory  optimization  schemes.  With  the  addition  of  recent  advancements 
in  the  field  of  trajectory  optimization,  results  suitable  for  controlling  a  full  6-DOF 
model  should  be  possible  without  needing  a  large  portion  of  those  simplifying  assump¬ 
tions. 

2.2  The  Trajectory  Optimization  Problem 

2.2.1  Problem  Formulation.  The  trajectory  optimization  problem,  also 
known  as  an  optimal  control  problem,  falls  under  the  broader  umbrella  of  Dynamic 
Optimization  and,  following  the  Hull’s  notation  format  [12],  is  characterized  by  the 
following  statement:  Determine  the  control  history  u(t )  that  minimizes  the  perfor¬ 
mance  index 


(2.1) 


subject  to  the  system  dynamics 


X  =  f(t,x,u), 


(2.2) 


the  specified  initial  conditions 


t(  0)  =  lo,  x(0)  =  x0> 


(2.3) 


the  specified  final  conditions 

V>( tf,xf )  =  0, 

(2.4) 

the  path  control  constraints 

C(t,  x,  u)  <  0, 

(2.5) 

and  the  path  state  constraints 

S(t,  x)  <  0. 

(2.6) 

The  performance  index,  J  in  Equation  2.1,  also  known  as  the  system’s  Cost 
Function,  is  a  scalar  function  of  the  system’s  state  vector,  x,  control  vector,  u,  time, 
t,  final  time,  tf,  and  state  vector  at  that  final  time,  Xf.  Since  this  research  focuses 
on  minimizing  the  time  required  to  complete  a  specified  maneuver,  the  performance 
index  can  be  simplified  to 

J  =  t,.  (2.7) 

Another  common  parameter  to  include  in  the  performance  index  would  be  the  fuel 
expended  over  the  course  of  the  maneuver.  This  would  create  a  minimum  fuel  opti¬ 
mization  setup. 

The  system  dynamics  in  Equation  2.2,  also  known  as  the  system’s  equations  of 
motion  (EOM),  are  those  equations  which  govern  how  the  aircraft  behaves  in  flight. 
As  depicted,  these  equations  are  also  functions  of  time,  the  current  state,  and  the 
control  inputs. 

The  initial  and  final  conditions  in  Equations  2.3  and  2.4  respectively  define 
the  state  variables  and  time  at  both  the  initiation  and  termination  of  the  maneuver. 
Unlike  the  dynamic  constraints,  where  all  variables  must  be  defined,  the  initial  and 
final  condition  functions  do  not  require  that  all  of  the  variables  be  defined,  just  enough 
to  fully  define  the  desired  maneuver.  Normally,  one  would  fully  define  the  initial  states 
and  time  and  only  define  those  final  states  and  time  that  are  necessary  to  define  the 
target  state.  For  example,  a  minimum  time  to  climb  problem,  where  you  were  not 
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concerned  with  any  of  the  final  state  variables  except  the  altitude  would  only  require 
a  definition  of  the  altitude  variable  for  the  final  conditions  constraints. 

The  control  and  state  inequality  path  constraints  included  in  equations  2.5  and 
2.6  define  the  operating  boundaries  of  the  system.  In  the  case  of  the  control  history, 
these  constraints  represent  the  physical  control  control  deflection  limits.  For  example, 
without  these  path  constraints,  the  optimization  system  may  try  to  push  a  throttle 
setting,  which  will  be  defined  later  as  0  <  5t  <  1,  beyond  its  limits  and  achieve  some 
unfeasible  result.  The  state  path  constraints  include  both  physical  limitations,  such 
as  defining  sea  level  as  a  minimum  altitude,  and  additional  factors  which  are  included 
to  define  the  exact  problem  that  the  user  is  trying  to  simulate.  By  coupling  the  ability 
to  change  the  desired  initial  and  final  conditions  with  the  ability  to  add  or  remove 
path  constraints  on  any  of  the  state  variables,  the  user  is  provided  a  robust  tool  which 
can  simulate  a  wide  variety  of  maneuvers. 

Each  of  these  parameters  is  combined  into  an  adjoint  cost  function  through  the 
use  of  Lagrange  Multipliers  which  are  also  known  as  the  adjoint  or  costate  variables  [3]. 
The  major  impact  of  this  is  that,  by  changing  any  of  the  the  basic  constraints  in  the 
problem,  the  cost  function  which  is  actually  being  optimized  also  changes.  As  a 
result,  the  solution  to  an  optimal  control  problem  is  highly  dependent  of  the  problem 
formulation  and  the  method  of  adjoining  the  constraints. 

2.2.2  Solution  Methods.  Given  this  basic  problem,  there  are  many  methods 
for  attempting  to  solve  the  trajectory  optimization  problem.  Many  of  these  meth¬ 
ods  are  covered  by  Betts  in  an  enlightening  survey  paper  [2],  In  this  paper,  Betts 
categorizes  the  majority  of  solution  methods  into  two  categories:  indirect  and  di¬ 
rect  methods.  Indirect  methods,  as  described  by  Betts,  require  the  determination 
of  explicit  solutions  of  the  equations  defined  by  the  problem’s  EOM,  the  necessary 
and  transversality  conditions,  and  the  maximum  principle.  In  effect,  these  methods 
don’t  solve  the  actual  optimal  control  problem,  but  instead  create  a  dualization,  or 
transformation,  of  that  problem  by  way  of  the  Hamiltonian  and  then  solve  that  new 
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problem  [22].  In  order  to  accomplish  this,  one  is  usually  required  to  solve  a  non-linear 
multipoint  boundary-value  problem.  If  this  description  weren’t  daunting  enough,  the 
utility  of  indirect  methods  is  fairly  limited  due  to  the  fact  that  one  must  not  only 
have  analytical  expressions  of  the  EOM,  but  also  possess  a  good  guess  as  to  the  fi¬ 
nal  answer.  Even  with  these  caveats,  indirect  methods  have  a  very  small  region  of 
convergence  [2],  For  these,  and  other,  reasons,  indirect  methods  are  rarely  used  for 
anything  other  than  fairly  simple  problems. 

With  direct  methods,  on  the  other  hand,  a  solution  is  found  by  manipulating 
parametric  representations  of  the  state  and/or  control  variables  to  directly  affect  the 
objective  function.  These  methods,  also  known  as  nonlinear  programming  problems 
(NLP),  are  generally  the  preferred  choice,  as  they  do  not  require  the  labor  intensive 
dualization  of  the  problem  and  the  analytic  derivation  of  the  necessary  conditions 
associated  with  that  task  [22],  Additionally,  direct  methods  generally  have  larger 
convergence  regions,  which  in  turn  creates  less  stringent  restrictions  on  the  initial 
guess.  Of  each  of  these  methods,  the  most  common  are  various  forms  of  what  are 
known  as  shooting  methods.  Shooting  methods  basically  take  initial  guesses  of  the 
optimal  control  histories  and  integrate  the  EOM  to  determine  the  performance  index 
associated  with  those  guesses.  As  expected,  these  methods  still  require  fairly  good 
initial  guesses  before  a  solution  is  even  feasible  and  even  then  can  result  in  massive 
numbers  of  iterations.  Unfortunately,  these  methods  generally  cannot  handle  what 
one  author  terms  “industry-strength”  problems.  These  problems,  as  is  often  the  case 
for  higher  end  aircraft  flight  simulations,  are  usually  characterized  by  complexities 
such  as  non-differentiable  table-lookups  [22], 

Over  the  last  several  years,  a  great  deal  of  work  has  gone  into  development  of 
advanced  trajectory  optimization  techniques.  This  push  has  resulted  from  increasing 
needs  in  two  major  areas:  satellite  orbit  transfers  and  Unmanned  Air  Vehicle  (UAV) 
control.  Work  in  the  area  of  UAV  control  has  ranged  from  efforts  looking  at  aircraft 
engagements  of  air  defense  systems  [15,16],  to  guidance  in  windy  environments  [14,28], 
to  real-time  trajectory  optimization  [18].  However,  in  each  of  these  cases,  in  order 
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to  find  solutions  using  traditional  trajectory  optimization  methods,  the  problems  are 
simplified  to  the  point  that  the  results,  though  often  very  enlightening,  can  not  be 
translated  into  useful  products  for  either  real  aircraft  applications  or  even  high  end  6- 
DOF  simulations.  Many  of  these  issues  can  now  be  addressed  due  to  recent  advances 
in  the  field  of  Pseudospectral  (PS)  methods. 

2.2.3  Legendre  Pseudospectral  Method.  Depending  on  which  terms  are  dis- 
critized,  direct  solution  methods  can  be  divided  into  several  categories,  the  most 
common  of  which  are  control  parameterization  and  state  and  control  parameteriza¬ 
tion  [13].  The  previously  mentioned  shooting  methods  are  examples  of  control  pa¬ 
rameterization  solution  methods  in  that  the  control  history  is  approximated  and  the 
differential  equations  are  propagated  forward  through  numerical  integration  schemes. 
As  would  be  expected,  the  size  of  the  NLP  significantly  increases  when  the  states  are 
included  in  the  parameterization,  but  state  and  control  parameterization  methods  are 
able  to  avoid  several  of  the  pitfalls  of  the  common  direct  shooting  methods. 

Pseudospectral  methods,  also  known  as  orthogonal  collocation  methods,  are  a 
subclass  of  state  and  control  parameterization  methods  which  approximate  the  states 
and  controls  with  a  finite  set  of  interpolating  polynomials  [13].  These  polynomials, 
evaluated  at  N  discritization  points  (nodes),  are  then  differentiated  and  constrained 
to  equal  the  differential  equations  of  the  original  problem  thereby  approximating  the 
state  derivatives. 

The  PS  solution  software  package  used  for  this  research  is  based  on  recent  de¬ 
velopments  in  the  area  of  Legendre  Pseudospectral  Methods  (LPM),  a  further  subset 
of  basic  PS  methods.  Legendre  Pseudospectral  Methods  are  characterized  by  their 
approximation  of  the  system’s  states  by  a  basis  of  N  Lagrange  interpolating  poly¬ 
nomials,  Ci(i  =  1, . . .  ,N),  and  by  N  nodes,  the  placement  of  which  are  defined  by 
Legendre-Gauss-Lobatto  (LGL)  points  [13].  Using  this  new  format,  the  state  vector 
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is  approximated  as 


N 

x[t)  «  X(t)  =  Y  £i(r)X(Ti),  (2.8) 

i= 1 

and  the  state  derivatives  at  the  kth  node  are  subsequently  approximated  as 

N 

x(rk)  «  X(rfc)  =  Y  Ci(rk)X(Ti),  {k  =  l,...,K).  (2.9) 

2—1 

The  control  vector  is  similarly  approximated  as 

N 

u(r)  &  Ujr)  =  Ci{r)U(Ti).  (2.10) 

i= 1 

The  resulting  problem  is  then  solved  via  a  nonlinear  programming  method  based  on 
sequential  quadratic  programming  [13]. 

Since  the  dualization  of  a  problem  via  indirect  methods  and  the  discritization  of 
the  same  problem  via  direct  methods  are  not  necessarily  commutative  operations,  the 
major  benefit  of  the  LPM  is  that  it  is  one  of  only  two  methods  which  have  been  shown 
to  preserve  the  order  of  the  original  problem  [22],  The  implication  of  this  statement, 
which  is  a  very  watered  down  form  of  the  Covector  Mapping  Principle,  is  that  the 
solution  to  an  LPM  problem  also  satisfies  the  problem’s  necessary  conditions,  thereby 
making  the  solution  method  both  direct  and  indirect  at  the  same  time  [22], 

2.2.4  Revised  Problem  Formulation.  In  order  for  a  PS  method  to  work, 
the  system’s  dynamics  and  other  equations  which  were  originally  used  to  define  the 
problem  must  be  translated  into  algebraic  constraint  equations  which  can  then  be 
applied  at  each  of  the  nodes  along  the  approximated  trajectory.  For  this  reason,  it 
is  desirable  to  rewrite  the  problem  formulation  in  terms  of  constraint  equations.  As 
described  by  Ross  [20],  the  problem  statement  is  now  the  following:  Determine  the 
state  and  control  pair,  {x,u\,  and  possibly  event  times,  To  and  77,  that  minimize  the 
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performance  index 


rTs 

J  =  E(x0,Xf,TQ,Tf)  +  /  F(t,x,  u)dr, 


subject  to  the  dynamic  constraints 


(2.11) 


x  =  }{t,x,u),  (2.12) 

the  event  constraints 

eL  <  e(x0,  xe,  xf,  t0,  re,  rf)  <  eu ,  (2.13) 

the  state  and  control  path  constraints 

hL  <  h(x,u,r )  <  hu ,  (2.14) 

and  the  state  and  control  variable  box  constraints 

xL  <  x(t)  <  xu ,  (2.15) 

uL  <  u(t )  <  uu ,  (2.16) 

The  cost  function,  J,  in  equation  2.11  is  known  as  a  Bolza  function  and  is 
comprised  of  an  event  cost  function,  E,  called  the  Mayer  Cost,  and  a  running  or 
integral  cost,  F,  called  the  Lagrange  Cost.  Note  that  this  problem  formulation  has 
benefited  from  a  variable  substitution  in  that  time  has  been  effectively  removed  from 
the  system  and  can  now  be  included  as  another  variable  that  can  be  solved  for.  This 
is  a  crucial  step  in  enabling  the  minimum  time  problem  which  is  the  main  focus  of 
this  research. 

The  dynamic  constraints  are  equivalent  to  the  previously  defined  dynamic  equa¬ 
tions  with  the  addition  of  a  r  substitution.  The  event  constraints,  e,  are  used  to  define 
the  boundary  and  internal  node  constraints  are  defined.  A  basic  problem  will  have 
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two  event  constraint  equations,  one  for  the  initial  conditions,  eo,  and  one  for  the  ter¬ 
minating  conditions,  e/,  as  described  previously.  Note  that  the  superscript  L  and  U 
are  used  to  signify  the  lower  and  upper  constraint  boundaries.  It  is  also  convenient  to 
introduce  the  idea  of  a  knot  at  this  point.  Any  node  along  the  trajectory  which  will 
have  specific  event  constraints  applied  to  it  will  be  referred  to  as  a  knot.  While  the 
initial  and  final  nodes  are  obviously  knots  at  time  r0  and  Tf  respectively,  a  knot  can 
also  be  defined  at  any  internal  event  time,  r0  <  re  <  Tf  as  well.  This  method  of  speci¬ 
fying  internal  knots,  incidently  known  as  the  knotting  method  [23],  is  extremely  useful 
in  defining  points  which  a  trajectory  must  pass  through  or  points  where  dynamics 
change,  such  as  a  multi-stage  rocket  losing  a  stage. 

The  path  constraints,  as  previously  defined  by  Equations  2.5  and  2.6,  are  now 
divided  into  two  constraint  categories.  This  division  is  performed  to  increase  numer¬ 
ical  efficiency  as  well  as  readability.  Those  path  constraints  which  are  imposed  over 
the  entire  trajectory,  state  and  control  physical  limits  for  example,  are  now  termed 
Box  Constraints  as  in  Equations  2.15  and  2.16.  Neither  path  nor  box  constraints  are 
required  to  fully  define  a  problem,  but  box  constraints  allow  a  user  to  confine  the 
solution  space  to  realistic  values,  provide  arbitrary  numerical  limits  to  ease  computa¬ 
tion  times,  and  avoid  known  singularities  in  the  mathematical  model.  For  example, 
instead  of  defining  the  Northing  Position,  Pjv,  limits  as  —  oo  <  P/v  <  oo,  which  is 
perfectly  reasonable,  one  would  define  the  limits  with  values  which  are  much  smaller 
than  infinity,  but  larger  than  any  value  which  the  system  is  expected  to  see.  In  this 
way,  the  constraints  do  not  actually  inhibit  a  solution,  but  do  drastically  decrease 
the  solution  space  over  which  the  optimization  routine  must  search.  This  method  of 
formulating  the  problem  through  a  series  of  constraints  will  be  readdressed  in  Chapter 
3  when  the  problem  for  this  specific  research  project  is  defined. 

2.2.5  Scaling  and  Balancing.  A  crucial  step  in  any  numerical  method  is 
to  properly  scale  and  balance  a  problem.  Without  properly  scaling  and  balancing 
a  problem  it  is  entirely  possible  that  an  otherwise  well  defined  problem  can  behave 
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very  poorly  within  a  numerical  method  routine.  The  basic  idea  behind  this  step  is  to 
convert  the  units  of  a  problem  from  a  physically  meaningful  unit  system  to  a  system 
which  behaves  in  a  desirable  fashion  during  an  optimization  routine  [11].  There  are 
two  basic  rules  for  transforming  the  unit  system  of  an  optimization  problem. 

First,  each  of  the  parameters  to  be  optimized  must  be  of  similar  magnitude  [11]. 
This  is  where  the  scaling  step  comes  into  play.  When  the  constraint  and  cost  equations 
of  an  optimization  problem  are  combined  into  an  adjoint  cost  function,  the  optimiza¬ 
tion  scheme  is  essentially  looking  at  scalar  representation  of  the  combination  of  each 
of  those  terms.  If  the  parameters  are  of  significantly  differing  orders  of  magnitude 
the  behaviors  of  one  parameter  might  mask  those  of  another,  possibly  more  impor¬ 
tant,  parameter.  For  a  simple  example  of  this  issue,  consider  a  satellite  in  an  orbit 
around  the  Earth  which  is  concerned  with  observing  some  property  of  the  Sun.  Two 
parameters  which  could  easily  be  included  in  this  problem  are  the  distance  between 
the  satellite  and  the  Sun  and  the  attitude  angles  of  the  satellite.  In  this  case,  one 
term  is  on  the  order  of  1  x  108  km  and  the  other  might  be  on  the  order  of  20°.  The 
most  common  solution  to  this  problem  is  to  scale  each  respective  parameter  by  a 
representative  value  in  a  linear  fashion  as  defined  by 


x  =  Dx  (2.17) 

where  x  is  scaled  version  of  the  state  vector,  x,  and  D  is  the  associated  scaling  matrix. 

After  the  problem’s  units  have  been  properly  scaled  there  still  exists  a  problem 
that  the  values  may  still  have  a  widely  varying  range.  Using  the  same  satellite  example 
from  before,  the  scaled  distance  from  the  satellite  to  the  Sun  might  vary  by  15,000  km , 
which  would  be  1.000015  in  the  scaled  unit  system,  but  an  attitude  angle  could  still 
vary  from  0°  to  360°,  which  would  be  0  to  18  in  a  dimensionless  systems  scaled  by  20°. 
The  relative  values  still  pose  a  problem  for  the  optimization  routine.  The  solution 
to  this  problem  is  to  balance  the  unit  system  by  shifting  the  units  such  that  each 
parameter  has  a  similar  range  [11].  A  common  method  for  performing  this  operation 


16 


is  to  use  the  minimum,  a,  and  maximum,  b,  values  which  a  parameter,  Xi,  is  expected 
to  see  during  the  course  of  the  trajectory  and  balance  the  basic  unit  according  to 


2  Xi  a,i  +  bi 


(2.18) 


which  can  also  be  written  as 

x  =  Dx  +  c.  (2-19) 

Through  the  use  of  Equation  2.19,  each  of  the  parameters  in  the  optimization  routine 
now  exist  in  the  same  range  as  defined  by  —  1  <  x$  <  1. 


2.3  Equations  of  Motion 

2.3.1  Summary  of  Assumptions.  Mathematical  models  of  real  systems  are, 
at  best,  approximations  of  the  actual  systems.  In  almost  all  situations,  it  is  im¬ 
practical,  and  often  impossible,  to  fully  capture  every  aspect  of  a  physical  system  in  a 
mathematical  representation.  As  discussed  in  Section  2.1,  the  requirement  for  making 
significant  simplifying  assumptions  before  a  problem  can  be  attempted  is  a  prevalent 
issue  in  aircraft  agility  research.  Though  this  effort  avoids  major  assumptions,  such 
as  the  assumption  of  a  point  mass,  there  are  several  minor  assumptions  which  were 
made  in  the  construction  of  the  mathematical  model.  The  assumptions  inherent  to 
the  mathematical  model  used  in  this  research  are  listed  in  Table  2.1.  The  assumptions 

Table  2.1:  Summary  of  the  assumptions  inherent  in  the  mathematical  model. 

1.  The  aircraft  is  assumed  to  be  a  rigid  body. 

2.  The  aircraft  is  assumed  to  be  symmetric  about  the  x-z  plane. 

3.  The  aircraft  does  not  expend  fuel  and  therefore  has  a  constant  mass. 

4.  Control  actuator  dynamics  are  neglected. 

5.  No  external  wind  or  environment  effects  are  included. 

6.  The  Earth  is  assumed  to  be  flat  and  non- rotating. 

7.  Gravity  is  assumed  to  be  constant  throughout  the  entire  reference  frame. 

in  Table  2.1  include  common  assumptions  about  the  physical  properties  of  the  aircraft 
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as  well  as  the  properties  of  the  physical  environment.  This  collection  of  assumptions, 
though  fairly  numerous,  does  not  alter  the  problem  significantly. 


Figure  2.1:  Aircraft  Axis  Definitions. 


2.3.2  Coordinate  systems.  The  development  of  the  aircraft  equations  of 
motion,  as  used  in  this  research,  required  three  basic  coordinate  systems.  The  first 
of  these  is  the  basic  aircraft  body-axis  system.  This  system,  as  seen  in  Figure  2.1,  is 
defined  as  having  it’s  origin  at  the  aircraft’s  center  of  gravity  (CG)  with  the  x-axis 
out  the  nose  of  the  aircraft,  the  y- axis  out  the  right  wing,  and  the  £-axis  out  the 
bottom  of  the  aircraft  as  defined  in  a  standard  right-handed  system. 

The  second  coordinate  system  used  is  the  stability-axis  system.  The  stability- 
axis  system,  as  also  seen  in  Figure  2.1,  is  a  rotation  of  the  body-axis  system  about 
the  y-axis,  by  an  angle  of  a,  where  a  is  the  angle  between  the  x-axis  in  the  body- axis 
system  and  the  relative  wind.  The  direction  cosine  matrix  (DCM)  associated  with 
this  rotation  is  defined  as  follows: 


cos  (a) 

0 

sin(a) 

C^  = 

0 

1 

0 

(2.20) 

—  sin(o:) 

0 

cos(a) 
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where  the  subscripts  b  and  s  are  respectively  used  to  denote  the  body-axis  and 
stability-axis  systems.  It  should  be  noted  that  the  main  reason  for  writing  the  equa¬ 
tions  of  motion  in  the  stability-axis  coordinate  system,  as  will  be  done  later,  is  that 
it  offers  some  fairly  significant  advantages  when  attempting  to  linearize  the  system. 
Though  linearization  is  not  used  in  this  research,  the  aircraft  model  used  was  previ¬ 
ously  developed  for  those  purposes  and  for  reasons  of  consistency,  the  equations  of 
motion  were  kept  in  the  stability-axis  coordinate  system. 

Finally,  an  inertial  reference  system  was  required  for  the  purpose  of  tracking 
the  relative  position  of  the  aircraft  throughout  a  maneuver  since  neither  the  body- 
axis  system  nor  the  stability-axis  system  is  suitable  for  this  purpose.  The  standard 
North- East- Down  (NED)  navigation  coordinate  system  was  used  for  this  purpose. 
This  is  an  earth  fixed  coordinate  system  with  the  origin  at  an  arbitrary  point  on  the 
surface  of  the  Earth  ((0,0,0)  in  the  case  of  this  research)  with  the  x-axis  aligned 
with  the  North  direction,  the  y-axis  aligned  with  the  East  direction,  and  the  2-axis 
pointed  into  the  Earth  as  defined  in  a  standard  right-handed  system.  In  this  case,  the 
system  is  also  defined  as  the  inertial  reference  system.  The  relationship  between  the 
navigation  system  and  the  body  system  is  defined  by  the  Euler  angles  in  a  standard 
yaw-pitch-roll  sequence,  as  seen  in  the  following  DCM: 


c9  op 

c9  sip 

—s9 

Cbn  = 

( — ccj)  sip  +  scp  s9  cip) 

( ccj)  cip  +  scp  s9  sip) 

scf)  c9 

(2.21) 

(scp  sip  +  ccp  s9  cip ) 

(—scp  cip  +  ccj)  s9  sip) 

ccp  c9 

where  the  subscript  n  is  used  to  denote  the  navigation  system.  Also  note  that  for 
space  saving  purposes  in  Equation  2.21  the  trigonometric  functions  sine  and  cosine 
are  abbreviated  as  s  and  c  respectively. 

As  would  be  expected  from  the  definition  of  the  body-axis  system,  Figure  2.2 
depicts  the  positive  directions  for  each  of  the  six  degrees  of  freedom  modeled  in  the 
equations  of  motion.  The  three  translational  degrees  of  freedom  are  denoted  as  posi- 
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Figure  2.2:  Degree  of  Freedoms  Definitions. 


tive  in  the  directions  of  the  positive  axes  and  the  three  rotational  degrees  of  freedom 
are  defined  as  positive  according  to  right  handed  rotations  about  the  respective  axes. 

2.3.3  Aircraft  Equations  of  Motion.  The  aircraft  model  used  for  this  effort 
is  based  on  the  F-16  model  presented  by  Stevens  and  Lewis  [26].  The  code  that 
was  used  to  simulate  this  model  is  contained  in  a  set  of  MATLAB  files  created  by 
previous  AFIT  students  as  replacements  for  the  FORTRAN  files  provided  by  Stevens 
and  Lewis.  The  original  FORTRAN  files  as  well  as  the  derived  MATLAB  files  have 
been  used  for  numerous  years  and  have  been  extensively  tested  to  ensure  correct 
results.  The  only  modifications  made  to  the  model  for  this  effort  are  all  found  in  the 
physical  properties  of  the  aircraft  that  will  be  described  in  Chapter  4.  The  baseline 
dynamics  and  equations  of  motion  remained  the  same. 

The  aircraft  model  is  defined  by  the  state  vector 
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\ 

VT 

/ 

Velocity 

a 

Angle  of  Attack 

Sideslip  Angle 

<f> 

Euler  Roll  Angle 

9 

Euler  Pitch  Angle 

i> 

Euler  Yaw  Angle 

Ps 

>  =  < 

Roll  Rate 

Qs 

Pitch  Rate 

Rs 

Yaw  Rate 

Pn 

Northing  Position 

Pe 

Easting  Position 

h 

Vertical  Position 

pow 

Engine  Thrust  Dynamics  Lag  State 

Similarly,  the  control  vector  is  defined  as: 


(2.22) 


\ 

5x 

Throttle  Control 

8e 

Elevator  Control 

>  =  < 

5a 

Aileron  Control 

5r 

Rudder  Command 

(2.23) 


It  should  be  noted  here  that  the  aircraft’s  altitude,  h,  is  defined  as  the  negative  of 
the  vertical  position  in  the  NED  system  and  that  the  aircraft’s  pitch  rate  (Qs)  in  the 
stability  system  is  exactly  equal  to  the  pitch  rate  (Q/,)  in  the  body-axis  system  since 
the  transformation  is  about  the  y- axis  and  therefore  does  not  affect  the  pitch  rate. 
It  is  only  denoted  as  Qs  here  for  purposes  of  consistency.  In  terms  of  the  controls, 
positive  control  deflections  are  defined  as  those  which  cause  negative  rotations  about 
their  respective  axes  as  previously  defined.  Using  this  convention,  the  equations  of 
motion  in  the  stability  axis  are  found  in  Table  2.2. 
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Tabic  2.2:  Aircraft  Equations  of  Motion  [26]. 


FORCE  EQUATIONS 


VT 


Ft  cos(a  +  ap  cos (/?)  —  D  +  rrigi 
m 


— Ft  sin  (a  +  aT)  —  L  +  1115:3  +  m Vr(Q  cos (/3)  —  Ps  sin(/3)) 

rnVU  cos(/3) 

•  — Ft  cos(a  +  ay)  sin  /3  —  C  +  11152  —  mVrPs 

111  Vt 


(2.24) 

(2.25) 

(2.26) 


KINEMATIC  EQUATIONS 

4>  =  P  +  tan(6l)  (Q  sin(0)  +  R  cos(</>) ) 

6  =  Q  cos(0)  —  Psin(0) 

•  Q  sin(0)  +  R  cos(0) 

'*’= - ^ - 


(2.27) 

(2.28) 

(2.29) 


MOMENT  EQUATIONS 


P 


Jxz[Jx  —  Jy  +  Jz\PQ  —  [Jz(Jz  —  Jy )  +  +  J zl  +  Jxn 


Jr.  J 7 


-  72 

°xz 


r,  (Jz  -  Jx)PR  -  Jxz(P2  —  R2)  +  m 

y  J 

Jy 

A  \{Jx  Jy)Jx  F  Jxz\PQ  J xz\J x  Jy  F  Jz\QR  F  Jxz^  F  Jx 77. 

*  ~  J  J  -  J2 

■'xz 

NAVIGATION  EQUATIONS 

Pat  =  UcO  cp  +  V (— c0  s-0  +  S0  s#  ci/j )  +  IU (s0  s-0  +  c0  c-0) 
Pe  =  Uc9  sp  +  V  ( c9  cijj  F  s</>  s9  sp)  +  W  (— scj)  cp  +  ccj)  sO  sp) 
h  =  Us9  —  Vs<j>  c9  -  Wccj)  c9 


(2.30) 

(2.31) 

(2.32) 

(2.33) 

(2.34) 

(2.35) 
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There  are  several  variables  and  notations  included  in  Table  2.2  for  space  and 
ease  of  reading  purposes  which  must  be  defined  for  the  reader.  In  the  Force  Equations, 
Equations  2.24-2.26,  L ,  D,  and  C  respectively  define  the  Lift,  Drag,  and  Side  forces 
acting  on  the  aircraft.  In  this  same  equation,  the  Thrust  Force  and  Thrust  Angle 
are  represented  by  Ft  and  ax  respectively.  Additionally,  Ps  is  the  roll  rate  in  the 
stability  axis  and  Rs  is  similarly  the  yaw  rate  in  the  stability  axis.  Finally,  the  wind 
axis  gravity  terms  are  defined  as: 

gi  =  g(—cac/3s9  +  s/3s(pc9  +  sac/3c(j)c9 ) 

()2  =  g(cas/3s9  +  c/3s<fic9  —  sas/3c(pc9 )  (2.36) 

#3  —  g(sas9  +  cac(j)c9 ) 

The  moment  equations,  Equations  2.30-2.32,  contain  the  moment  of  inertia 
( Jx , Jy ,  Jz)  and  cross-product  of  inertia  ( Jxz ,  Jxy,  Jyz)  terms.  Note  that  here  it  is 
assumed  that  the  x-z  plane  is  a  plane  of  symmetry  for  the  aircraft,  which  causes  each 
of  the  cross-products  of  inertia  terms,  except  Jxz,  to  be  equal  to  zero.  The  moment 
equations  also  include  the  torque  terms  Z,  m,  and  n.  It  is  through  these  torque  terms, 
as  depicted  in  Equations  2.37-2.39, 


l  =  f{P,P,R,6a,6r), 

(2.37) 

f  (  5  ^ ^  ?  Q  ?  &T  7  ^ e  5  Ta  ) , 

(2.38) 

=  f(P,P,R,6a,6r,Tf,), 

(2.39) 

that  the  elevator,  aileron,  and  rudder  are  able  to  exert  their  control  over  the  aircraft 
[26].  Note  that  Ta  and  Tp  in  Equations  2.38  and  2.39  are  the  thrust  offset  angles, 
both  of  which  are  assumed  to  be  zero  for  this  research.  Additionally,  the  Navigation 
Equations,  Equations  2.33-2.35,  contain  the  velocity  components  U,  V,  and  W  in  the 
body-axis  system.  These  components  of  velocity  are  related  to  Vr,  a,  and  (3  through 
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Equations  2.40-2.42. 


VT  =  VU2  +  V2  +  W2 


(2.40) 


a  =  tan  1  (2-41) 

/?  =  sin-1  (2.42) 

The  addition  of  the  engine  thrust  dynamics  lag  state,  pow,  is  an  artifact  of  the 
specific  aircraft  model  used  in  this  research.  As  will  be  discussed  later,  this  state  will 
be  represented  as  a  first-order  lag  as  in  Equation  2.43, 


pow  =  —  ( powc  —  pow), 


(2.43) 


where  T  is  the  lag  constant  and  powc  is  the  commanded  power  setting  which  is  a 
function  of  the  commanded  throttle  setting,  5t- 
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III.  Modeling  &;  Implementation 

3.1  Aircraft  Models 

Three  aircraft  models  were  used  for  this  endeavor,  all  of  which  are  based  on 
the  NASA-Langley  F-16  wind  tunnel  test  data  presented  by  Stevens  and  Lewis  [26] 
and  the  mathematical  model  derived  from  that  data.  In  order  to  demonstrate  the 
difference  in  time-optimal  trajectories  for  various  aircraft,  three  derivatives  of  the 
basic  F-16  model  were  utilized.  The  three  models  used  include  a  baseline  F-16  model, 
a  model  with  increased  thrust,  and  a  model  with  increased  wing  area.  These  changes 
were  chosen  as  a  way  of  modifying  the  two  most  important  aircraft  performance 
parameters:  thrust-to-weight  ratio  ( T/W )  and  wing  loading  ( W/S ). 

With  all  else  being  equal,  an  aircraft  with  higher  T/W  will  be  able  to  reach 
higher  velocities,  climb  faster,  accelerate  quicker  and  sustain  higher  turn  rates  [19]. 
On  the  other  hand,  an  aircraft  with  a  lower  W/S  will  be  capable  of  lower  stall  speeds 
and  tighter  instantaneous  turns  [19].  Together,  the  three  models  should  simulate  a 
baseline  fighter  aircraft,  a  second  fighter  which  is  more  prone  to  vertical  maneuvers 
and  a  third  aircraft  which  is  more  prone  to  horizontal  maneuvers. 

3.1.1  Physical  Layout.  The  F-16,  as  depicted  in  Figure  3.1,  is  a  classic  single 
engine  multi-role  tactical  aircraft  design.  The  design  features  a  single  vertical  tail  and 
a  wing  and  stabilizer  configuration  which  creates  a  neutrally  stable  platform.  Control 
of  the  aircraft  is  obtained  through  the  use  of  the  ailerons,  elevators,  and  rudder.  Note 
that  the  actual  F-16  elevators  are  differentially  controllable,  which  allows  them  to 
roll  the  aircraft,  especially  during  transonic  speeds.  The  mathematical  model  does 
not  include  the  ability  to  differentially  articulate  the  elevators  which  means  that 
the  elevators  cannot  be  used  for  roll  control.  The  aircraft  has  a  single  afterburning 
turbofan  engine  along  the  centerline  which  is  assumed  to  act  along  the  x-axis  in  the 
body  system.  The  configuration  in  Figure  3.1  is  the  nominal  configuration  of  the 
baseline  aircraft  model.  The  two  model  variants  developed  in  during  this  research 
effort  are  notional  only  and  are  therefore  not  physically  depicted  here. 
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Figure  3.1:  Aircraft  Layout  [1]. 

3.1.2  Physical  Parameters.  The  relevant  physical  parameters  for  the  base¬ 
line  aircraft  model  are  summarized  in  Table  3.1.  These  values  are  obtained  from  the 
mathematical  model  description  provided  by  Stevens  and  Lewis  in  Reference  [26]  and 
are  all  assumed  to  remain  constant  throughout  each  simulation.  The  only  value  in 
Table  3.1  which  deviates  from  the  Stevens  and  Lewis  model  is  the  value  of  the  CG 
location.  The  original  CG  location  value  was  actually  0.35  c,  but,  as  will  be  discussed 
in  Chapter  4,  it  was  necessary  to  move  the  CG  location  to  0.25  c  for  stability  reasons. 

Table  3.1:  Summary  of  Baseline  Aircraft  Model  Physical  Parameters 


Parameter  Value 


Weight  (W) 

20,500 

lbs 

Jx 

9,496 

slug  ' 

■ ft 2 

Jy 

55,814 

slug  ' 

■ft2 

JZ 

63,100 

slug  ' 

■ft2 

Jx  z 

982 

slug  ' 

■ft2 

CG  Location 

25 

%  MAC 

Span  (6) 

30 

ft 

Area  (S) 

300 

ft 2 

MAC  (c) 

11.32 

ft 
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As  previously  mentioned,  the  thrust-to-weight  ratio  and  wing  loading  of  the  two 
notional  models  were  varied  by  increasing  the  available  thrust  or  increasing  the  wing 
area  of  the  model.  Table  3.2  summarizes  the  changes  to  the  baseline  model  which 
were  used  to  create  the  two  notional  models.  Note  that  Max  Thrust  is  defined  as  the 
static  thrust  available  at  sea  level.  These  deviations  are  equivalent  to  a  25%  increase 
in  available  thrust  for  Model  //  2  and  a  25%  increase  in  wing  area  for  Model  //  3. 


Table  3.2:  Summary  of  Aircraft  Models  Deviation  from  Baseline  Model 


Parameter 

Model  #  1 

Model  #  2 

Model  #  3 

Weight  (W) 

20,500 

20,500 

20,500 

lbs 

Max  Thrust  (T) 

20,000 

25,000 

20,000 

lbs 

Wing  Area  (S) 

300 

300 

375 

ft2 

Wing  Loading  {W / S') 

68.33 

68.33 

54.67 

lbs/ ft2 

Thrust-to-Weight  Ratio  (' T/W ) 

0.9756 

1.2195 

0.9756 

- 

3.1.3  Aerodynamic  Model.  The  aerodynamic  model  consists  of  a  series  of 
table  lookup  routines  which  are  used  to  calculate  the  aerodynamic  force  and  moment 
buildup.  These  calculations  are  all  completed  in  the  body-axis  frame  and  then  trans¬ 
ferred  to  the  stability-axis  system  and  navigation  system  as  required.  A  summary  of 
the  aerodynamic  lookup  tables  is  found  in  Table  3.3.  The  sign  convention  for  each  of 
the  coefficients  adhere  to  the  conventions  outline  in  Figure  2.2.  A  separate  routine  is 
also  provided  to  model  the  standard  atmosphere. 

The  data  in  these  lookup  tables  was  modified  from  the  original  NASA-Langley 
data  by  Stevens  and  Lewis  to  include  the  effects  of  the  F-16’s  leading-edge  flap.  The 
data  for  the  leading-edge  flap,  the  actual  deployment  of  which  is  scheduled  based 
on  Angle-of- Attack  (AOA)  and  Mach  Number,  was  originally  contained  in  several 
other  lookup  tables.  By  merging  the  data  tables,  Stevens  and  Lewis  were  able  to 
significantly  reduce  the  number  of  table  lookups  required  while  still  maintaining  the 
leading-edge  flap  effects  minus  the  associated  actuator  dynamics  [26]. 
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Table  3.3:  Summary  of  Aerodynamic  Lookup  Tables. 


Coefficient 


Component  Function  of 


Basic 

Basic 

Basic 


a 

oi ,  5e 
(3,  5a,  $r 


Damping  Derivatives 
X-axis  Force  (Cx) 
Y-axis  Force  (Cy) 

Z-axis  Force  (Cz) 
Rolling  Moment  (Ci) 

-  Aileron  Component 

-  Rudder  Component 
Pitching  Moment  (Cm) 
Yawing  Moment  (Cn) 

-  Aileron  Component 

-  Rudder  Component 


Basic 

,  5a ,  5, 

Basic 

a,  j3 

5a 

a,  /3 

5r 

a,  (3 

Basic 

a,  5e 

Basic 

a,  [3 

5a 

a,  [3 

5r 

a,  f3 

3.1.4  Engine  Model.  The  engine  model  is  based  on  a  model  developed 
by  NASA-Langley  in  parallel  with  their  wind  tunnel  measurement  efforts.  Again, 
this  data  is  provided  by  Stevens  and  Lewis  [26].  The  engine  model  is  comprised  of 
two  MATLAB  routines.  The  first  routine  is  a  model  of  the  engine  power  response 
to  the  power  setting  commanded  and  is  modeled  as  a  basic  first-order  lag  with  a 
variable  time  constant  which  is  a  function  of  the  actual  engine  power  setting  (pow) 
and  the  commanded  power  setting  ( powc ).  The  second  routine  is  a  table  lookup  of 
the  thrust  values  as  a  function  of  the  power  setting  ( pow ),  the  aircraft’s  altitude  (h), 
and  the  current  Mach  Number  (M).  The  lookup  tables  include  data  for  idle  power, 
military  power,  and  maximum  power,  where  the  change  from  military  power  setting 
to  maximum  power  settings  occur  when  pow  >  0.77  [26].  Also  of  note  is  that  the 
engine  is  modeled  as  having  a  constant  angular  momentum  of  160  slug  ■  ft2 / s. 

3.1.5  Model  Control  Actuator  Limits.  Within  the  mathematical  model,  the 
control  surface  actuators  are  modeled  with  both  deflection  and  rate  limits.  Addition¬ 
ally,  they  are  modeled  with  a  basic  first-order  lag.  A  summary  of  the  limits  on  the 
control  system  is  provided  in  Table  3.4. 


Table  3.4:  Summary  of  control  limits  inherent  in  mathematical  model  [26]. 

Control  Deflection  Limit  Rate  Limit  Time  Constant 

Elevator  ±25.0°  60°/s  0.0495  s  lag 

Aileron  ±21.5°  80°/s  0.0495  s  lag 

Rudder  ±30.0°  120°/s  0.0495  s  lag 

3.1.6  Model  State  Limits.  As  will  be  shown  later,  the  solution  space  for 
the  trajectory  optimization  problem  is  defined  by  the  physical  limits  on  state  and 
control  variables.  Before  the  optimization  routine  will  work  correctly,  there  are  two 
key  sets  of  limits  which  must  be  addressed:  those  which  are  defined  by  the  physics  of 
the  problem,  and  those  which  are  inherent  in  the  mathematical  model.  Limits  which 
are  defined  by  the  physics  of  the  problem  are  fairly  intuitive.  For  example,  except  for 
a  few  rare  situations,  none  of  which  were  considered  here,  the  altitude  of  an  aircraft 
is  limited  by  h  >  0. 

On  the  other  hand,  limits  inherent  in  the  mathematical  model,  though  arguably 
more  important  to  the  optimization  routine,  are  not  nearly  as  obvious  without  a  good 
working  knowledge  of  the  model  being  used.  These  limits  correspond  to  the  limits 
of  the  model  itself  and  are  the  result  of  various  assumptions  which  were  made  in 
the  construction  of  the  mathematical  model.  In  the  case  of  the  F-16  model  used  in 
this  effort,  the  majority  of  these  limits  are  the  result  of  incomplete  data  sets  within 
table  lookups.  For  example,  without  an  intimate  knowledge  of  the  aircraft  model, 
one  would  not  know  that  the  altitude  is  actually  restricted  to  0  ft  <  h  <  50,  000  ft 
by  the  simple  fact  that  these  are  the  limits  of  the  available  altitude  data  from  the 
original  wind  tunnel  testing. 

Table  3.5  contains  a  summary  of  the  limits  on  the  state  variables  which  are 
inherent  in  the  F-16  mathematical  model.  One  item  of  note  within  this  table  is  the 
upper  limit  on  AO  A.  The  data  lookup  tables  support  angles  as  high  as  45°,  but  the 
upper  limit  was  lowered  to  25°  to  avoid  stall  and  the  deep  stall  issues  inherent  in 
the  F-16  model  [26].  Allowing  the  optimization  routine  to  search  in  this  region  would 
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simply  waste  time  and  possibly  confuse  the  algorithm.  The  remaining  limits  are  based 
on  the  limits  of  the  data  lookup  tables  for  various  components  within  the  model. 


Table  3.5:  Summary  of  state  limits  inherent  in  mathematical  model. 


Parameter 


Minimum  Maximum  Units 


Mach  Number  (M) 
Altitude  (h) 


0 

0 

-10 

-30 


1 

50,000  ft 

45  degrees 

30  degrees 


Angle  of  Attack  (a) 
Sideslip  Angle  (f3) 


3.2  Optimization  Software 

Solutions  for  the  complex  trajectory  optimization  problems  in  this  research  were 
possible  through  the  use  of  a  new  and  novel  dynamic  optimization  software  package, 
which  was  created  by  researchers  at  the  Naval  Post-Graduate  School  and  published  by 
Elissar,  LLC.  This  package,  which  is  provided  as  a  set  of  MATLAB  p-code,  exploits 
the  advantages  of  Pseudospectral  Methods  to  create  a  user  friendly  and  versatile 
optimization  package. 

The  DIDO  software  package,  named  after  Queen  Dido  of  Carthage  who  was 
the  first  person  to  solve  a  dynamic  optimization  problem  [21],  is  a  robust  dynamic 
optimization  package  which  is  specifically  tailored  to  the  optimal  control  problem. 
Through  the  use  of  this  package,  a  user  is  able  to  define  an  optimal  control  problem 
in  a  very  intuitive  fashion. 

3.2.1  Problem  Definition  Structure.  Within  the  DIDO  framework,  trajec¬ 
tory  optimization  problems  are  defined  by  sets  of  equality  and  inequality  constraints. 
These  constraints  fully  define  the  system  dynamics,  the  valid  solution  space,  and  the 
performance  index  for  the  optimization  problem.  This  setup  is  accomplished  through 
the  use  of  five  specific  hies,  which  are  basically  the  only  input  to  the  optimization 
routine. 
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3.2. 1.1  Problem  Setup  File.  The  first,  and  most  important,  of  the 
DIDO  input  hies  is  the  Problem  Setup  File.  This  hie  is  used  to  define  the  numerical 
values  of  the  scaling  factors  and  constraints  for  each  of  the  state  and  control  variables. 
In  this  capacity,  the  values  in  this  hie  completely  define  a  specihc  maneuver  and 
this  hie  is  therefore  the  only  hie  which  needs  to  be  changed  to  simulate  a  different 
maneuver.  To  this  end,  the  setup  hie  contains  the  initial  and  target  state  values  as 
well  as  any  state  values  which  are  prescribed  along  the  trajectory  path  as  described 
in  Equation  2.13. 

Additionally,  the  setup  hie  contains  the  limit  values  for  the  path  and  box  con¬ 
straints  from  Equations  2.14-2.16.  Unlike  the  other  constraint  equations,  which  are 
explicitly  defined  in  the  additional  input  hies,  the  box  constraint  equations  are  as¬ 
sumed  to  be  in  the  form  of  Equations  2.15  and  2.16  and  DIDO  enforces  these  con¬ 
straints  internally. 

It  should  be  noted  that  this  structure  is  a  slight  deviation  from  the  format  sug¬ 
gested  in  the  DIDO  User’s  Manual  [21],  in  which  the  numerical  constraint  information 
is  all  included  in  the  problem’s  Main  hie.  This  change  was  made  due  to  the  simple  fact 
that,  with  a  model  comprised  of  13  states  and  4  controls,  the  amount  of  information 
needed  to  define  the  problem  warranted  its  own  hie  for  readability  purposes. 

3.2. 1.2  Dynamic  Constraints  File.  The  second  DIDO  input  hie,  the 
dynamic  constraints  hie,  contains  the  problem’s  differential  equations  and  is  very 
similar  in  structure  to  the  input  function  for  the  various  initial  value  problem  solvers, 
such  as  ode23 ,  within  MATLAB.  The  dynamics  hie  receives  the  values  of  the  state  and 
control  variables  at  a  specihc  point  along  a  notional  trajectory  and  then  calculates 
the  differential  values  of  the  state  variables  and  returns  those  values  to  DIDO.  DIDO 
then  compares  this  information  to  what  it  predicted  the  differentials  to  be  based  on 
the  derivatives  of  the  interpolating  polynomials  and  then  updates  that  prediction 
accordingly. 
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In  this  capacity,  the  dynamic  constraints  hie  is  required  to  contain  all  of  the 
functions  and  equations  necessary  to  determining  the  derivatives  of  the  state  variables. 
For  this  research,  this  means  that  the  dynamic  constraints  hie  is  the  equivalent  of  the 
F  hie  presented  in  Stevens  and  Lewis  [26]  and  contains  the  entire  mathematical  model 
of  the  aircraft  motion  dynamics  and  all  of  the  requisite  table  lookup  calls. 

One  consideration  for  this  hie  is  that  the  dynamics  hie  must  contain,  or  be 
able  to  retrieve,  the  system’s  scaling  factors.  This  is  necessary  for  situations  where 
the  dynamic  equations  are  written  in  the  unsealed  unit  system  since  the  input  states 
provided  by  DIDO  are  already  scaled.  Similarly,  the  state  derivatives  returned  to 
DIDO  must  also  be  returned  to  the  scaled  unit  system. 

3.2. 1.3  Event  Constraints  File.  The  Events  hie  contains  the  event 
constraint  functions  laid  out  in  Equation  2.13  as  functions  of  the  state  variables  along 
the  trajectory.  Note  that  the  Events  hie  does  not  contain  the  actual  target  values 
for  the  desired  state  points  as  previously  laid  out  in  the  Problem  Setup  hie,  just  the 
equations  needed  to  determine  how  closely  the  trajectory  meets  those  requirements. 

In  the  case  of  maneuvers  which  are  dehned  by  an  initial  state  and  a  hnal  target 
state,  the  Events  hie  contains  a  function  for  each  of  the  state  variables  which  are 
prescribed  at  the  initial  state  and  another  set  of  functions  for  each  of  the  state  variables 
that  are  used  to  dehne  the  target  state.  Any  additional  points  which  the  trajectory 
must  pass  through  would  also  be  dehned  in  the  Events  hie. 

3. 2. 1-4  Cost  Constraints  File.  The  Cost  Constraint  hie  is  where  the 
user  defines  the  cost  function  which  will  be  used  for  the  problem.  This  simple  hie 
breaks  the  Bolza  cost  function  from  Equation  2.11  into  the  sum  of  the  single  point 
Mayer  cost  function  and  the  integral  Lagrange  cost,  thus  allowing  the  user  to  dehne 
an  event  cost  as  well  as  a  running  cost. 

3. 2. 1.5  Path  Constraints  File.  Similar  in  structure  to  the  Events  hie, 
the  Path  Constraints  hie  contains  the  path  constraint  equations  from  Equation  2.14. 
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This  file  is  an  optional  addition  to  the  problem  since  it  is  not  always  necessary  to 
restrict  the  path  of  the  states  or  controls. 

3.2.2  Method  Verification.  In  order  to  verify  both  that  the  DIDO  software 
worked  as  advertised  and  that  it  was  implemented  correctly,  several  problems  with 
known  solutions  were  simulated  for  comparison  purposes.  The  most  representative  of 
these  comparisons  was  a  sample  problem  posed  and  solved  by  Bryson  [6].  The  problem 
(Problem  4.5.24  [6]),  is  to  find  the  Minimum  Time  to  Climb  for  a  727  aircraft  climbing 
2000  ft  from  sea  level  and  returning  to  its  initial  velocity  and  flight  path  angle.  The 
equations  of  motion  for  the  aircraft,  in  normalized  units,  are  simplified  to  the  following 
four  equations: 


T  cos(a  +  e)  —  Cr>V2  —  sin (7), 

(3.1) 

:  T  sin(a  +  e)  —  ClV2  —  cos(7), 

(3.2) 

h  =  V  sin(7), 

(3.3) 

x  —  V  cos(7), 

(3.4) 

where  V,  7,  and  h  are  the  state  variables  and  a  is  the  only  control  variable  in  this 
free  final  time  problem. 

The  results,  as  presented  by  Bryson,  are  shown  in  Figure  3.2(a)  and  (b)  where 
the  solid  blue  lines  are  the  actual  response  and  the  dashed  red  lines  are  the  steady- 
state  optimal  solution.  These  results,  which  match  those  provided  in  the  book,  were 
found  by  using  the  sample  MATLAB  script  hies  provided  by  Bryson  and  the  internal 
MATLAB  function  fmincon.  On  the  other  hand,  solving  the  same  problem  through 
the  use  of  DIDO  provides  the  results  found  in  Figure  3.2(c)  and  (d). 

Comparing  the  results  from  the  baseline  Bryson  solutions  and  the  DIDO  output 
reveals  a  few  interesting  discrepancies.  First,  though  the  results  are  very  similar,  the 
optimum  time  values  provided  by  the  two  methods  do  not  match.  As  seen  in  Table  3.6, 
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(a)  Bryson  State  and  Control  Trajectory 


(b)  Bryson  Aircraft  Energy  Response 


(c)  DIDO  State  and  Control  Trajectory  (d)  DIDO  Aircraft  Energy  Response 


(e)  Propagated  State  and  Control  Trajectory 


(f)  Propagated  Aircraft  Energy  Response 


Figure  3.2:  Results  for  the  Bryson  727  sample  optimization  problem. 


the  DIDO  solution  results  in  a  final  time  over  a  full  second  less  than  that  predicted 
by  Bryson  through  the  use  of  fmincon. 
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Parameter  fmincon  DIDO 

Ff  56.114  55.003  sec~ 

Table  3.6:  Bryson  727  problem  solution  comparison. 

Since,  as  previously  discussed,  the  results  from  DIDO  are  an  approximation  of 
the  actual  solution,  it  is  possible  that  the  results  are  a  poor  approximation  of  the 
actual  aircraft  trajectory.  The  simple  way  to  check  that  the  approximation  is  valid  is 
to  propagate  the  aircraft  states  over  time  using  a  numerical  integration  scheme.  The 
results  of  propagating  the  states  forward  using  MATLAB’s  odef5  function  are  found 
in  Figure  3.2(e)  and  (f).  The  propagated  results  match  the  approximate  trajectories 
almost  exactly,  which  means  that  the  aircraft  is  actually  capable  of  reaching  the  target 
state  by  following  the  trajectory  predicted  by  DIDO. 

The  next  obvious  question  is:  why  are  the  DIDO  results  better  than  the  fmincon 
results?  The  differences  between  the  results  are  caused  by  a  conglomeration  of  several 
different  factors.  First,  as  noted  by  Bryson  in  the  comments  in  his  sample  code  [6], 
this  problem  has  very  poor  global  convergence  properties  and  converges  very  slowly 
when  not  near  what  it  deems  to  be  the  optimal  solution.  Second,  the  guess  for 
the  DIDO  problem  setup  and  the  fmincon  setup  are  different.  The  guess  given  in 
the  fmincon  is  actually  the  optimal  control  history  and  the  final  solution  is  limited 
to  be  within  10%  of  the  guess.  This  is  not  a  good  way  to  determine  an  optimal 
solution,  but  this  was  the  way  the  code  was  provided  by  Bryson,  most  likely  as  the 
result  of  numerous  iterations,  and  deviations  from  this  method  result  in  the  algorithm 
converging  to  a  different  local  minimum.  The  guess  of  the  control  history  provided  to 
DIDO  is  a  constant  a  =  0,  which  is  neither  realistic  nor  representative  of  the  actual 
solution.  This  was  done  to  demonstrate  the  robustness  inherent  within  DIDO  since 
giving  a  similar  guess  in  the  fmincon  method  results  not  only  in  it  failing  to  find  the 
optimal  solution,  but  in  it  actually  failing  to  converge  to  any  solution.  Additionally, 
the  solution  found  from  using  DIDO  is  also  replicated  if  the  initial  guess  is  reverted 
back  to  the  guess  provided  by  Bryson.  This  shows  that  the  pseudospectral  method  is 
not  as  susceptible  to  falling  into  non-optimal  local  minimums  as  the  fmincon  method. 
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3.3  Simulation  Setup 


3.3.1  Scaling  and  Balancing.  As  previously  mentioned,  the  choice  of  scaling 
and  balancing  factors  can  significantly  affect  the  results  of  an  optimization  routine. 
To  that  end,  each  state  and  control  parameter  was  converted  to  a  unitless  system 
with  typical  values  ranging  from  —  1  <  x  <  1  through  the  use  of  Equation  2.19.  The 
scaling  and  balancing  factors  chosen  to  accomplish  this  are  listed  in  Table  3.7.  These 
terms  were  chosen  based  on  the  range  of  values  which  each  parameter  is  expected  to 
see  for  a  typical  maneuver.  Note  that  the  boundaries  of  the  typical  values  are  not 
necessarily  the  same  as  the  box  constraints,  which  will  be  addressed  later. 


Table  3.7:  Scaling  and  balancing  factors  used  in  the  optimization  scheme. 


Parameter 

Typical  Values 

Lower  Upper 

Scaling  Factor 

Balancing  Factor 

(Unit  less) 

V 

0 

1000 

ft/s 

500 

-1 

a 

-10° 

25° 

deg 

17.5° 

-0.4286 

P 

-15° 

15° 

deg 

15° 

0 

1 

00 

o 

o 

180° 

deg 

180° 

0 

6 

-87° 

87° 

deg 

87° 

0 

■0 

-180° 

180° 

deg 

00 

o 

o 

0 

P 

I 

00 

o 

o 

180° 

deg/ s 

o 

o 

00 
T— I 

0 

Q 

-45° 

45° 

deg/ s 

45° 

0 

R 

-45° 

45° 

deg/ s 

45° 

0 

Pn 

-25000 

25000 

ft 

25000 

0 

Pe 

-25000 

25000 

ft 

25000 

0 

h 

10000 

30000 

ft 

10000 

-2 

pow 

0 

100 

50 

-1 

5t 

0 

1 

0.5 

-1 

-24° 

24° 

deg 

24° 

0 

da 

-21.5° 

21.5° 

deg 

21.5° 

0 

Sr 

-30° 

30° 

deg 

30° 

0 

t 

0 

100 

s 

50 

-1 

3.3.2  Box  Constraints.  One  of  the  most  critical  steps  in  solving  the  trajec¬ 
tory  optimization  problem  for  a  fairly  complicated  6-DOF  mathematical  model,  was 
the  development  of  the  box  constraints.  As  previously  discussed,  the  box  constraints 
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must  include  the  limitations  of  the  model,  the  physical  limits,  and  reasonably  large 
bounds  on  otherwise  unbounded  parameters.  Table  3.8  lists  the  box  constraints  uti¬ 
lized  in  the  Problem  Setup  File  in  both  the  original  and  scaled  and  balanced  unit 
systems. 


Table  3.8:  Box  constraints  used  in  the  optimization  scheme. 


Parameter 

Standard  Units 

Lower  Upper 

Designer  Units 

Lower  Upper 

V 

0.1 

1000 

ft/s 

-0.9998 

1 

a 

-10° 

25° 

deg 

-1 

1 

P 

o 

O 

CO 

1 

30° 

deg 

-2 

2 

0 

-540° 

Oi 

o 

o 

deg 

-3 

3 

9 

o 

IT- 

OO 

I 

o 

00 

deg 

-1 

1 

0 

-540° 

Oi 

o 

o 

deg 

-3 

3 

P 

-270° 

270° 

deg  /  s 

-1.5 

1.5 

Q 

1 

00 

o 

o 

00 

o 

o 

deg/ s 

-4 

4 

R 

o 

o 

00 
T— 1 

I 

00 

o 

o 

deg/ s 

-4 

4 

Pn 

-100,000 

100,000 

ft 

-4 

4 

Pe 

-100,000 

100,000 

ft 

-4 

4 

h 

0 

50,000 

ft 

-2 

3 

pow 

0 

100 

-1 

1 

St 

0 

1 

-1 

1 

Se 

1 

to 

4^ 

o 

24° 

deg 

-1 

1 

Sa 

-21.5° 

21.5° 

deg 

-1 

1 

Sr 

o 

O 

CO 

1 

30° 

deg 

-1 

1 

t 

0.001 

250 

s 

-0.99998 

4 

There  are  several  important  pieces  of  information  which  fed  into  the  develop¬ 
ment  of  these  constraints.  First,  and  foremost,  were  the  model  limitations  imposed 
by  the  table  lookup  data  in  the  aerodynamic  and  engine  models.  These  limits  defined 
the  box  constraints  for  V,  a,  (3,  h,  pow ,  St,  Se,  Sa  and  Sr.  The  lower  limits  of  V  and 
t  are  chosen  as  small  numbers  which  were  relatively  close  to  zero  without  creating 
the  singularity  which  an  actual  value  of  zero  would  cause.  Similarly,  the  constraints 
on  6  were  chosen  to  avoid  the  inherent  singularity  at  ±90°.  The  remaining  limits  on 
0,  0,  P,  Q,  R,  Pn,  Pe,  and  the  upper  limit  for  t  were  chosen,  after  experimentation 
with  the  model,  as  reasonably  large  values  which  are  used  in  place  of  ±oo  in  order 
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to  place  reasonable  limits  on  the  solution  space  and  subsequently  reduce  the  required 
run  time  for  the  optimization  scheme  without  impacting  the  final  solution. 

3.3.3  Maneuvers.  In  order  to  evaluate  the  effectiveness  and  utility  of  the 
trajectory  optimization  approach  in  determining  aircraft  agility  characteristics,  a  wide 
range  of  maneuvers  were  simulated.  The  simulated  maneuvers  have  been  divided  into 
three  main  categories:  Demonstration  Maneuvers,  Agility  Maneuvers,  and  Compound 
Maneuvers.  Each  of  the  maneuvers  starts  from  an  initial  state  with  the  aircraft 
trimmed  for  steady,  wings  level  flight  with  the  spatial  setup  defined  in  Table  3.9.  The 
actual  trimmed  values  of  the  remaining  states  and  controls  are  found  through  the  use 
of  a  trimmer  routine  modified  from  Stevens  and  Lewis  [26]  and  vary  slightly  with  each 
aircraft  due  to  their  thrust  and  surface  area  differences. 

Table  3.9:  Initial  spatial  setup  for  each  maneuver. 

Parameter  Initial  Condition 
V  500  Jtfs 

PN  0  ft 

Pe  0  ft 

h  20,000  ft 

From  this  initial  state,  each  of  the  three  aircraft  models  is  tasked  to  perform 
each  maneuver  in  a  time  optimal  fashion  for  later  comparison  and  analysis. 

3.3.3. 1  Demonstration  Maneuvers.  The  maneuvers  in  the  Demon¬ 
stration  Maneuvers  category  are  used  as  an  initial  demonstration  of  the  utility  of 
trajectory  optimization.  Through  the  three  maneuvers  in  this  category,  the  basic 
validity  of  the  approach  is  demonstrated  and  verified  before  moving  on  to  the  more 
rigorous  maneuvers  in  the  Agility  Maneuvers  category. 

Northing  Position  Change 

The  first  maneuver  in  the  Demonstration  category  is  the  Northing  Position 
Change.  In  this  maneuver,  each  aircraft  is  tasked  to  move  downrange  by  10,000  ft  as 
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fast  as  possible.  The  basic  setup  is  intended  to  resemble  the  classic  Brachistochrone 
problem  since  this  is  a  very  common  initial  example  in  optimization  literature.  The 
brachistochrone  problem,  as  posed  by  John  Bernoulli  in  1696  [7],  is  to  find  the  shape 
of  a  frictionless  wire  which  would  allow  a  bead  sliding  along  it  to  move  from  one  point 
to  another  while  only  being  acted  upon  by  gravity.  The  solution  is  a  shape  known  as 
a  cycloid,  which  is  defined  as  a  path  which  is  generated  by  a  point  on  a  circle  that  is 
rolling  in  a  horizontal  direction  without  slipping  [7]. 

Though  the  Northing  Position  Change  Maneuver  differs  from  the  classic  brachis¬ 
tochrone  maneuver  through  the  inclusion  of  the  aircraft  dynamics  as  well  as  an  initial 
velocity,  the  maneuver  is  intentionally  designed  to  produce  a  similar  trajectory.  The 
maneuver  is  numerically  defined  by  two  events:  the  initial  state,  and  the  target  state. 
Starting  from  the  initial  state  discussed  previously,  the  aircraft  is  tasked  to  transition 
to  the  state  described  in  Table  3.10.  Note  that  the  variables  with  a  target  value  of 
”  are  considered  to  be  free  variables  and  are  allowed  to  take  on  any  value  within 
the  solution  space  at  that  point.  For  this  maneuver,  the  only  restriction  on  the  final 
state  is  that  the  aircraft  must  be  10,000  ft  North  of  its  original  position. 

Table  3.10:  Numerical  definition  of  the  Northing  Position  Change  maneuver. 

State  Variable  Initial  State  Target  State 


V 

500 

- 

ft/s 

a 

C^trim 

- 

deg 

P 

0 

- 

deg 

0 

- 

deg 

e 

@trim 

- 

deg 

if 

0 

- 

deg 

p 

0 

- 

deg/ s 

Q 

0 

- 

deg/ s 

R 

0 

- 

deg/ s 

Pn 

0 

10,000 

ft 

Pe 

0 

- 

ft 

h 

20,000 

- 

ft 

pow 

pOW^rim 

- 
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Unconstrained  3-D  Position  Change 

The  Unconstrained  3-D  position  Change  maneuver,  as  defined  in  Table  3.11,  is 
intended  to  demonstrate  an  out-of-plane  maneuver.  For  this  maneuver,  each  aircraft 
is  tasked  to  transition  to  a  point  10,000  ft  North  and  East  of  its  original  position  and 
return  to  the  original  altitude.  The  remaining  states  are  left  as  free  variables. 

Table  3.11:  Numerical  definition  of  the  Unconstrained  3-D  Position  Change  ma¬ 
neuver. 

State  Variable  Initial  State  Target  State 


V 

500 

- 

ft/s 

a 

C^trim 

- 

deg 

(3 

0 

- 

deg 

<t> 

0 

- 

deg 

e 

@trim 

- 

deg 

fp 

0 

- 

deg 

p 

0 

- 

deg/ s 

Q 

0 

- 

deg/ s 

R 

0 

- 

deg / s 

Pn 

0 

10,000 

ft 

Pe 

0 

10,000 

ft 

h 

pow 

20,000 

POW  trim 

20,000 

ft 

Constrained  3-D  Position  Change 

Building  on  the  previous  maneuver,  the  Constrained  3-D  Position  Change  ma¬ 
neuver  tasks  the  aircraft  with  performing  the  same  basic  maneuver,  but  further  con¬ 
strains  the  target  state.  As  defined  in  Table  3.12,  this  maneuver  requires  that  the 
aircraft  move  to  a  point  10,000  ft  North  and  East  of  its  initial  position  and,  aside 
from  the  position,  return  to  its  original  wings  level,  trimmed  state. 

3. 3. 3. 2  Agility  maneuvers.  After  demonstrating  the  utility  and  ro¬ 
bustness  of  the  trajectory  optimization  scheme,  the  maneuvers  in  the  Agility  Maneu¬ 
vers  category  are  intended  to  be  more  rigorous  tests  of  the  limits  of  each  aircraft’s 
capabilities. 
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Table  3.12:  Numerical  definition  of  the  Constrained  3-D  Position  Change  maneuver. 


State  Variable  Initial  State  Target  State 


V 

500 

V0 

ft/s 

a 

^trinn 

«o 

deg 

P 

0 

A) 

deg 

0 

0 

00 

deg 

e 

@trim 

0o 

deg 

0 

0 

00 

deg 

P 

0 

P0 

deg  /  s 

Q 

0 

Qo 

deg  /  s 

R 

0 

Ro 

deg/ s 

Pn 

0 

10,000 

ft 

Pe 

0 

10,000 

ft 

h 

20,000 

ho 

ft 

pow 

pOWtrim 

pow0 

Bank  Angle  Capture 

The  first  agility  maneuver,  the  Bank  Angle  Capture  maneuver,  is  one  which  is 
widely  suggested  as  a  possible  agility  metric.  This  maneuver,  as  defined  in  Table  3.13, 
requires  that  the  aircraft  achieve  and  hold  a  90°  bank  angle.  To  make  the  maneuver 
more  rigorous,  the  aircraft  is  also  required  to  return  to  its  initial  heading,  velocity, 
and  cross-range  position. 


Unconstrained  Heading  Capture 

Similar  to  the  Bank  Angle  Capture  maneuver,  the  Unconstrained  Heading  Cap¬ 
ture  maneuver  is  intended  to  test  the  nose  pointing  capabilities  of  an  aircraft.  For 
this  maneuver  the  aircraft  is  required  to  reorient  its  nose  to  a  direction  defined  by 
6  =  0°  and  0  =  45°.  Though  the  aircraft  must  capture  this  attitude,  there  are  not 
any  requirements  for  the  aircraft  to  actually  be  moving  in  that  specific  direction.  In 
that  context,  this  maneuver  is  intended  to  simulate  a  situation  where  a  pilot  would 
like  to  rapidly  reorient  the  nose  of  the  aircraft  to  a  point  where  a  weapon  may  be 
employed  as  quickly  as  possible.  The  numerical  definition  of  this  maneuver  is  found 
in  Table  3.14. 
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Table  3.13:  Numerical  definition  of  the  Bank  Angle  Capture  maneuver. 


State  Variable  Initial  State  Target  State 


V 

500 

V0 

ft/s 

a 

C^trim 

- 

deg 

f> 

0 

- 

deg 

0 

0 

90° 

deg 

e 

@trim 

- 

deg 

0 

0 

00 

deg 

P 

0 

P0 

deg  /  s 

Q 

0 

Qo 

deg/ s 

R 

0 

Ro 

deg/ s 

Pn 

0 

- 

ft 

Pe 

0 

Pe0 

ft 

h 

20,000 

- 

ft 

pow 

pOWfrim 

- 

Table  3.14:  Numerical  definition  of  the  Unconstrained  Heading  Capture  maneuver. 

State  Variable  Initial  State  Target  State 


V 

500 

- 

ft/s 

a 

Ottrim 

- 

deg 

P 

0 

- 

deg 

0 

0 

- 

deg 

e 

@trim 

0 

deg 

0 

0 

45° 

deg 

P 

0 

P0 

deg/ s 

Q 

0 

Qo 

deg/ s 

R 

0 

Ro 

deg/ s 

Pn 

0 

- 

ft 

Pe 

0 

- 

ft 

h 

20,000 

- 

ft 

pow 

pOWtrim 

- 

Constrained  Heading  Capture 

Taking  the  Unconstrained  Heading  Capture  a  step  further,  the  Constrained 
Heading  Capture  maneuver  requires  that  the  aircraft  change  its  course  to  a  heading 
of  -0  =  90°  and  return  to  its  original  steady  level  flight  conditions.  The  numerical 
definition  for  the  Constrained  Heading  Capture  maneuver  is  found  in  Table  3.15. 
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Table  3.15:  Numerical  definition  of  the  Constrained  Heading  Capture  maneuver. 


State  Variable  Initial  State  Target  State 


V 

500 

V0 

ft/s 

a 

C^trim 

«o 

deg 

P 

0 

A) 

deg 

0 

0 

00 

deg 

e 

@trim 

0o 

deg 

ip 

0 

O 

O 

o 

deg 

p 

0 

Po 

deg  /  s 

Q 

0 

Qo 

deg/ s 

R 

0 

Ro 

deg  /  s 

Pn 

0 

- 

ft 

Pe 

0 

- 

ft 

h 

20,000 

ho 

ft 

pow 

pOWtrim 

pow0 

Position-Free  Heading  Reversal 

The  Position- Free  Heading  Reversal  is  the  first  of  three  heading  reversal  maneu¬ 
vers  designed  to  fully  tax  the  capabilities  of  each  aircraft.  Similar  to  the  Bank  Angle 
Capture  maneuver,  the  generic  heading  reversal  maneuver  is  often  suggested  as  a  air¬ 
craft  agility  metric.  Additionally,  previous  investigation  of  this  specific  maneuver  by 
Bocvarov  [4]  provides  a  good  means  of  comparing  the  full  6-DOF  results  with  results 
found  through  the  use  of  a  point-mass  simplifying  assumption. 

In  this  maneuver,  as  defined  in  Table  3.16,  the  aircraft  is  tasked  to  reverse  its 
heading  and  recapture  the  initial  steady  level  flight  conditions  in  minimum  time.  To 
that  end,  the  only  free  variables  in  this  maneuver  are  the  three  position  variables  Pn, 
Pe  and  h. 


Position-Fixed  Heading  Reversal 

The  Position-Fixed  Heading  Reversal  further  constrains  the  target  state  of  the 
Position-Free  Heading  Reversal  maneuver.  In  this  maneuver,  the  aircraft  is  required 
to  fully  recapture  its  initial  state  with  the  only  variation  being  the  new  heading  angle 
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Table  3.16:  Numerical  definition  of  the  Position-Free  Heading  Reversal  maneuver. 


State  Variable  Initial  State  Target  State 


V 

500 

V0 

ft/s 

a 

C^trim 

«o 

deg 

P 

0 

A) 

deg 

0 

0 

00 

deg 

e 

@trim 

0o 

deg 

ip 

0 

180° 

deg 

p 

0 

Po 

deg  /  s 

Q 

0 

Qo 

deg/ s 

R 

0 

Ro 

deg  /  s 

Pn 

0 

- 

ft 

Pe 

0 

- 

ft 

h 

20,000 

- 

ft 

pow 

pOWtrim 

pow0 

of  ip  —  180°.  For  comparison  purposes,  this  maneuver  was  also  previously  investigated 
by  Bocvarov  [4], 

Table  3.17:  Numerical  definition  of  the  Position-Fixed  Heading  Reversal  maneuver. 

State  Variable  Initial  State  Target  State 


V 

500 

V0 

ft/s 

a 

^ trim 

«o 

deg 

P 

0 

Po 

deg 

0 

0 

00 

deg 

9 

@trim 

0o 

deg 

ip 

0 

180° 

deg 

P 

0 

Po 

deg  /  s 

Q 

0 

Qo 

deg  /  s 

R 

0 

Ro 

deg  /  s 

Pn 

0 

Pn0 

ft 

Pe 

0 

Pe0 

ft 

h 

20,000 

ho 

ft 

pow 

POWtrim 

pow0 

Position-Free  Heading  Reversal  with  Altitude  Floor 

A  third  twist  on  the  Heading  Reversal  maneuvers  is  the  inclusion  of  a  restrictive 
minimum  altitude.  As  will  be  shown  later,  the  inclusion  of  an  altitude  floor  will  only 
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have  a  large  impact  on  the  Position-Free  Heading  reversal  maneuver,  so  this  is  the  only 
one  which  is  fully  investigated.  The  numerical  definition  for  this  maneuver  is  exactly 
the  same  as  that  found  in  Table  3.16  for  the  Position-Free  Heading  Reversal.  The 
difference  between  the  two  maneuvers  is  that  the  Box  Constraints  for  this  maneuver 
have  been  altered  such  that  the  minimum  altitude  is  raised  to  20,000  ft,  which  means 
that  the  aircraft  is  not  allowed  to  loose  any  altitude  during  this  maneuver. 


Initial  State  Capture 

Unlike  the  Position-Fixed  Heading  Reversal,  the  Initial  State  Capture  maneuver 
requires  the  aircraft  to  fully  recapture  its  initial  state  without  any  further  variations. 
This  means  that  the  aircraft  must  return  to  trimmed,  steady-level  flight  at  its  initial 
position,  attitude,  and  velocity.  The  numerical  definition  of  this  maneuver  is  found 
in  Table  3.18. 

Table  3.18:  Numerical  definition  of  the  Initial  State  Capture  maneuver. 

State  Variable  Initial  State  Target  State 


V 

500 

V0 

ft/s 

a 

O^trim 

«o 

deg 

P 

0 

A) 

deg 

0 

0 

00 

deg 

e 

@trim 

00 

deg 

ip 

0 

00 

deg 

p 

0 

Po 

deg/ s 

Q 

0 

Qo 

deg/ s 

R 

0 

Ro 

deg/ s 

Pn 

0 

Pn0 

ft 

Pe 

0 

Pro 

ft 

h 

20,000 

ho 

ft 

pow 

POW  trim 

pow  o 

Initial  State  Capture  with  Altitude  Floor 

As  will  be  shown  later,  the  Initial  State  Capture  maneuver  will  result  in  a 
trajectory  which  utilizes  an  altitude  loss  to  gain  energy  early  in  the  maneuver.  The 
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final  agility  maneuver,  the  Initial  State  Capture  with  Altitude  Floor,  is  designed 
to  remove  the  altitude  loss  option  from  the  Initial  State  Capture  maneuver.  The 
numerical  definition  for  this  maneuver  is  exactly  the  same  as  that  found  in  Table  3.18 
for  the  Initial  State  Capture.  The  only  difference  between  the  two  maneuvers  is  that 
the  Box  Constraints  for  this  maneuver  have  been  altered  such  that  the  minimum 
altitude  is  raised  to  20,000  ft ,  which  means  that  the  aircraft  is  not  allowed  to  loose 
any  altitude  during  this  maneuver. 

3. 3. 3. 3  Compound  Maneuvers.  Each  of  the  maneuvers  in  the  Demon¬ 
stration  and  Agility  categories  have  one  major  characteristic  in  common;  they  only 
look  at  the  transition  from  an  initial  state  to  a  final  state.  The  next  step  in  devel¬ 
oping  more  complicated  maneuvers  for  both  agility  evaluation  and  aircraft  control  in 
general  is  to  specify  a  series  of  states  which  the  aircraft  must  transition  through  on 
its  way  to  a  target  state.  With  this  observation,  a  discussion  is  warranted  on  the 
difference  between  a  trajectory  which  is  built  by  minimizing  the  time  to  each  point 
and  a  trajectory  which  is  built  by  minimizing  the  time  through  each  point. 

As  noted  by  Aides  [18],  the  minimum  time  trajectory  through  a  series  of  points 
is  almost  always  not  the  collection  of  minimum  time  trajectories  between  the  points 
in  question.  The  illustration  in  Figure  3.3  depicts  this  phenomenon  by  showing  two 
notional  minimum  time  trajectories  through  four  points.  In  the  first  case,  the  or¬ 
ange  line  represents  the  path  an  aircraft  might  follow  if  it  were  tasked  to  reach  each 
successive  point  in  minimum  time.  Traversing  from  the  initial  point  to  Point  A  in 
minimum  time  would  produce  a  straight  line  acceleration.  Moving  from  Point  A  to 
Point  B  in  minimum  time  would  ideally  be  a  straight  line  between  the  two,  but  the 
initial  trajectory  has  set  the  aircraft  up  with  a  large  velocity  in  the  direction  of  Point 
A  which  causes  the  aircraft  to  execute  a  turn  to  reacquire  a  path  to  Point  B.  This 
same  phenomenon  is  then  repeated  for  the  leg  from  Point  B  to  Point  C. 

On  the  other  hand,  if  the  aircraft  were  tasked  to  reach  Point  C  in  minimum  time 
while  passing  through  Points  A  and  B,  the  trajectory  might  look  something  like  the 
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maroon  line.  In  this  case,  the  aircraft  reaches  Point  A  after  the  first  aircraft,  but  by 
executing  a  heading  change  the  aircraft  has  set  itself  up  for  a  straight  line  acceleration 
through  Points  A  and  B  on  its  way  to  Point  C. 


—  Min  time  to  A,  B,  C 

Min  time  to  C,  through  A,  B 

Figure  3.3:  Notional  Minimum  Time  Trajectory 

4-Point  Position  Change 

To  demonstrate  the  multi-point  capabilities  inherent  in  DIDO  through  the  use 
of  the  knotting  method,  the  notional  maneuver  described  in  Figure  3.3  was  illustrated 
through  the  use  of  the  4- Point  Position  Change  maneuver.  As  defined  in  Table  3.19, 
the  aircraft  is  tasked  to  reach  a  final  state  in  minimum  time  while  passing  through 
two  separate  states  along  the  way.  The  setup  is  actually  identical  to  the  course  in 
Figure  3.3  in  that  the  final  three  points  all  lie  along  the  same  line  with  only  the 
starting  point  being  offset  from  that  line. 

3.3.4  Result  Verification.  The  final  step  in  determining  the  optimal  tra¬ 
jectories  is  to  verify  that  they  are  actually  feasible  trajectories.  Since  the  results 
returned  from  DIDO  are  only  approximations  of  the  state  and  control  trajectories, 
it  is  entirely  possible  the  results  will  be  poor  representations  of  the  actual  trajectory 
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Table  3.19:  Numerical  definition  of  the  4- Point  Position  Change  maneuver. 


State  Variable  Initial  State  Target  States 


A 

B 

Final 

V 

500 

- 

- 

- 

ft/s 

a 

Qtrim 

- 

- 

- 

deg 

P 

0 

- 

- 

- 

deg 

0 

- 

- 

- 

deg 

e 

@trim 

- 

- 

- 

deg 

ip 

0 

- 

- 

- 

deg 

p 

0 

- 

- 

- 

deg/ s 

Q 

0 

- 

- 

- 

deg/ s 

R 

0 

- 

- 

- 

deg/ s 

Pn 

0 

5,000 

5,000 

5,000 

ft 

Pe 

0 

0 

5,000 

10,000 

ft 

h 

20,000 

20,000 

20,000 

20,000 

ft 

pow 

POWtfifYi 

- 

- 

- 

of  the  aircraft.  In  order  to  verify  that  the  trajectories  are  achievable,  the  control 
trajectories  developed  by  DIDO  are  subsequently  used  as  the  input  to  a  numerical 
integration  scheme,  namely  ode45  in  MATLAB,  which  then  propagates  the  reactions 
of  the  aircraft  states  to  the  control  inputs.  The  resulting  state  trajectories  should 
closely  match  those  predicted  by  DIDO,  if  not,  then  the  results  are  not  valid. 
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IV.  Simulations  &;  Results 


4-1  Stability  Adjustments 

Initial  results  from  attempting  the  optimization  of  several  of  the  basic  maneuvers 
revealed  some  interesting  phenomena.  Initial  attempts  to  optimize  various  maneuvers 
not  only  took  an  exceedingly  long  amount  of  time  to  complete,  but  it  was  also  often 
the  case  that  the  optimization  results  would  not  match  the  propagated  results  at  all. 
Further  investigation  into  this  issue  revealed  two  instabilities  in  the  system:  one  in 
the  aircraft  model  and  another  in  the  optimization  scheme  itself. 

4-1.1  Aircraft  Model  Stability.  The  mathematical  model  of  the  F-16  used 
in  this  research  places  the  aircraft’s  center  of  gravity  at  the  aircraft’s  aerodynamic 
center  which  in  turn  creates  a  neutrally  stable  aircraft.  Although  this  is  an  accurate 
representation  of  the  actual  F-16,  a  real-life  F-16  utilizes  a  stability  augmentation 
system  to  keep  the  aircraft  from  going  unstable.  Since  the  basic  mathematical  model 
does  not  include  a  stability  augmentation  system,  even  fairly  small  control  deflections 
are  enough  to  cause  the  aircraft  to  go  unstable  and  depart  controlled  flight. 

In  general,  when  attempting  a  trajectory  optimization  problem,  the  issue  of  an 
inherently  unstable  aircraft  can  be  dealt  with  by  increasing  the  number  of  nodes  until 
they  occur  faster  than  the  frequency  of  the  instability.  Practically,  this  is  not  desirable 
due  to  the  fact  that  increasing  the  number  of  nodes  in  this  fashion  also  exponentially 
increases  the  time  required  for  convergence.  The  results  of  a  trade  study  revealed 
that  moving  the  aircraft’s  center  of  gravity  from  its  original  location  of  0.35c  to  a 
new  location  at  0.25c  would  artificially  stabilize  the  aircraft  model  without  greatly 
affecting  the  aircraft’s  maneuverability. 

4-1.2  Optimization  Routine  Stability.  The  second  instability  was  not  nearly 
as  easy  to  diagnose.  The  basic  symptom  of  this  problem  was  the  fact  that  the  state 
trajectories  resulting  from  the  propagation  of  the  aircraft  controls  did  not  match  the 
estimated  trajectories.  By  simply  looking  at  the  final  results  of  an  optimization  run, 
it  was  not  evident  what  was  causing  the  problem;  but  it  was  glaringly  evident  that  the 
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optimization  solution  was  incorrect.  This  basic  problem  is  illustrated  in  Figure  4.1 
where  the  blue  line  is  the  optimal  solution  provided  by  DIDO  and  the  red  line  is  the 
result  of  propagating  the  states  in  time  with  odeJ^5. 


Figure  4.1:  Departure  of  propagated  trajectory  from  the  optimal  solution. 

4- 1.2.1  Unstable  Responses.  Further  investigation  revealed  an  inter¬ 
esting  phenomenon  in  which  the  optimization  routine  actually  causes  the  unstable 
behavior.  This  phenomenon  is  visible  in  the  plots  in  Figure  4.2.  This  figure  contains 
a  control  history  and  a  3-D  trajectory  plot  for  each  of  three  points  in  the  optimization 
routine.  Though  defined  by  the  number  of  iterations  which  have  occurred  prior  to 
these  results,  the  specific  number  of  iterations  at  each  point  is  not  important  as  this 
phenomenon  occurs  at  different  points  in  the  iterative  process  for  different  maneuvers. 

Figure  4.2(a)  and  (b)  show  the  control  history  and  3-D  trajectory  for  the  Uncon¬ 
strained  3-D  Position  Change  maneuver  at  a  point  where  the  optimization  routine  has 
converged  to  a  fairly  good  solution.  In  the  trajectory  plot,  the  blue  line  represents  the 
approximate  trajectory  which  DIDO  thinks  it  is  following  and  the  red  line  represents 
the  “truth”  solution  obtained  by  propagating  the  controls  history  with  MATLAB’s 
ode45  numerical  integration  function.  In  this  particular  case,  the  red  and  blue  lines 
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lie  on  top  of  each  other,  which  means  that  the  approximation  is  fairly  good  at  this 
point. 

After  converging  to  a  decent  solution  at  30,000  iterations,  the  optimization 
routine  goes  to  work  on  whittling  down  the  cost  function,  which,  in  this  case,  is  the 
final  time.  By  100,000  iterations  the  routine  has  reduced  the  final  time  from  a  value  of 
22.838  s  at  30,000  iterations  to  a  new  value  of  21.885  s.  At  this  point  the  optimization 
routine  starts  looking  for  ways  to  shave  increasingly  small  amounts  of  time  off  of  the 
final  time. 

Figure  4.2(c)  depicts  the  result  of  the  routine’s  searching  for  ways  to  shave  off 
extra  time.  What  is  occurring  here  is  that  the  system  has  noticed  small  oscillations  in 
the  roll,  pitch,  and  yaw  response  of  the  aircraft  which  are  the  result  of  the  light  damp¬ 
ing  inherent  in  the  F-16  model.  Though  not  depicted  here  for  space  saving  purposes, 
there  are  many  examples  of  this  phenomenon  in  the  various  figures  in  Appendix  A. 
Basically,  the  system  begins  trying  to  nullify  those  small  oscillations  in  an  attempt  to 
smooth  the  trajectory  response.  In  doing  so,  the  system  begins  a  cycle,  similar  to  a 
Pilot  Induced  Oscillation,  where  it  begins  to  actually  drive  the  oscillations  and  then 
begins  to  increase  the  amplitude  of  the  controls  to  compensate  for  the  new  response. 

In  some  cases  the  system  is  able  to  recover  from  this  death  spiral,  but  most 
cases  result  in  a  situation  similar  to  Figure  4.2(e)  and  (f)  where  the  system  has  driven 
the  routine  unstable  and  caused  the  aircraft  to  depart  from  the  desired  trajectory. 
The  main  problem  here  is  that  the  optimization  system  does  not  know  that  the 
“optimal”  trajectory  is  no  longer  feasible  and  will  continue  along  this  path,  often 
actually  producing  results  which  it  thinks  are  optimal  solutions. 

1.2.2  Artificial  Stabilization  Methods.  Several  methods  of  artificially 
stabilizing  the  optimization  routine  were  attempted  with  varying  success.  By  modi¬ 
fying  the  problem’s  cost  function,  it  is  possible  to  effectively  penalize  the  routine  for 
control  histories  such  as  those  in  Figure  4.2.  When  this  phenomenon  was  first  discov¬ 
ered,  the  penalty  term  found  in  Equation  4.1  was  included  in  the  running  Lagrange 
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Aileron  (deg)  Throttle  (pet)  Aileron  (deg)  Throttle  (pet)  Aileron  (deg)  Throttle  (pet) 


(a)  Controls  at  30k  Iterations 


(b)  Trajectory  at  30k  Iterations 
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5  10  15  20 

Time  (s) 


(c)  Controls  at  100k  Iterations 


(d)  Trajectory  at  100k  Iterations 


(e)  Controls  at  200k  Iterations 


(f)  Trajectory  at  200k  Iterations 


Figure  4.2:  Optimization  induced  oscillation  development  over  time. 
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cost  term  (see  Equation  2.11)  of  the  performance  index. 


Pi  =  i(5i|5e|  +  B2\5a\  +  B3\Sr\),  (4.1) 

6tf 

The  basic  goal  of  this  term  is  to  penalize  the  system  for  the  large  magnitude  controls 
which  are  inherent  in  the  unstable  responses.  In  this  equation,  the  hat  over  the 
terms,  such  as  Se,  is  used  to  denote  that  these  are  values  in  the  scaled  and  balanced 
unit  system.  By  averaging  these  terms  and  dividing  by  the  final  time,  tf,  this  term 
is  reduced  to  a  scalar  value  where  0  >  P±  >  1.  This  scaling  makes  it  easier  to 
subsequently  scale  this  term  in  relation  to  the  Mayer  Cost  term  from  Equation  2.11 
to  adjust  the  relative  importance  of  the  final  time  and  the  penalty  term  through  the 
use  of  the  scaling  terms  Bn. 

By  itself,  the  inclusion  of  Equation  4.1  into  the  cost  function  resulted  in  faster 
convergence  for  a  majority  of  the  maneuvers  simulated,  but  still  couldn’t  handle 
certain  maneuvers  and  is  definitely  working  against  the  intent  of  pushing  the  limits 
of  the  aircraft’s  capabilities.  To  that  end,  a  second  Lagrange  parameter,  found  in 
Equation  4.2,  was  developed  to  address  the  issue  of  rapid  oscillations  in  the  controls 
instead  of  their  use  in  general. 

P2  =  -^(BiA5t  +  B2A8e  +  B3A8a  +  B4A8r),  (4.2) 

4  tf 

The  overall  intent  of  this  cost  function  is  to  force  the  system  to  smooth  the  control 
histories,  therefore  penalizing  the  cost  function  for  large,  rapid  changes  in  the  control 
history.  Once  again,  this  parameter  works  well  for  some  maneuvers  and  not  as  well 
for  others. 

The  next  iteration  of  the  penalty  term  borrows  a  page  from  the  Digital  Signals 
Processing  world  and  adapts  steady  state  statistical  methods  to  the  current  problem. 
In  essence,  the  new  cost  function  attempts  to  minimize  the  variance  (a2)  of  the  control 
signals,  which  again  is  an  attempt  to  smooth  the  control  histories.  The  problem  is 
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that  normal  methods  of  calculating  the  variance  of  a  signal  are  based  on  a  constant, 
or  stationary,  reference  signal  with  noise  [25] .  Since  the  control  histories  are  anything 
but  stationary,  a  slightly  different  approach  is  required.  By  breaking  the  signal  into 
small  pieces  and  analyzing  the  variance  of  that  section,  a  better  estimation  of  the 
variance  of  the  entire  signal  can  be  obtained  [25].  As  shown  in  Equation  4.3,  the  new 
cost  function  includes  the  average  of  the  variance  of  N  specific  sections  of  the  control 
histories. 

i  N 

P3  =  Atf  +  —  +  B2all(6e)  +  B3<j2N(8a)  +  B4a^r(Sr)),  (4.3) 

o 

Though  the  inclusion  of  this  cost  function  into  the  performance  index  stabilized  the 
system  for  every  maneuver  under  investigation  it  also  had  one  major  side  effect.  In 
order  to  increase  the  accuracy  of  the  approximation  of  the  signal’s  variance  it  is 
necessary  to  increase  the  number  of  sections  into  which  the  signal  is  divided.  The 
most  effective  means  of  accomplishing  this  was  to  create  a  “window”  of  interest  which 
would  be  propagated  along  the  control  signal  and  the  variance  of  the  signal  within 
that  “window”  would  be  calculated  at  each  step.  This  step  drastically  increased 
the  computational  time  required  to  complete  each  iteration  and  subsequently  each 
optimization  run. 

The  final  iteration  on  the  penalty  term  takes  the  basic  premise  of  Equation  4.3 
one  step  further  and  provides  the  basis  for  the  final  performance  index  used  through¬ 
out  the  remainder  of  this  research.  First,  a  noise  corrupted  stationary  signal  is  created 
by  taking  the  difference  between  a  control  signal  and  a  moving  average  MA  of  that 
signal.  The  variance  of  this  stationary  signal  is  used  as  the  variance  of  the  actual 
control  history  and  is  summed  with  the  variance  of  the  other  three  control  inputs. 
These  calculations  result  in  the  penalty  term  found  in  Equation  4.4. 

4 

Pi  =  £  -  A  (4.4) 

i= 1 


54 


The  inclusion  of  this  penalty  term  not  only  solved  the  instability  issues,  but  also 
drastically  reduced  the  computational  time  required  from  that  of  the  previous  cost 
function.  The  cost  function  used  throughout  the  remaining  portion  of  this  research  is 
seen  in  Equation  4.5.  The  relative  weighting  between  the  final  time  and  the  penalty 
terms  depict  the  relative  importance  of  each  term  to  the  final  solution. 

4 

J  =  loot/  +  10  (4.5) 

i= 1 

Unfortunately,  this  equation  also  depicts  the  fact  that  the  trajectory  optimization 
problem  which  uses  this  performance  index  is  no  longer  solving  for  a  time-optimal 
trajectory.  This  means  that  there  will  almost  always  exist  a  solution  which  would 
result  in  a  slightly  better  final  time.  Fortunately,  for  the  maneuvers  investigated  in 
this  research,  those  difference  are  negligible. 

The  cost  function  in  Equation  4.5  was  used  to  create  the  control  and  trajectory 
development  figures  in  Figure  4.3,  which  show  the  progression  of  the  controls  and 
trajectories  of  the  same  aircraft  and  maneuver  as  previously  demonstrated.  Note 
that  this  system  also  starts  to  exhibit  some  of  the  same  unstable  phenomenon  at  the 
100,000  iteration  mark.  The  difference  now  is  that  this  type  of  control  use  is  now 
penalized  and  the  system  converges  to  the  solution  in  Figure  4.3  (e)  and  (f). 

4-2  Results  Format 

The  results  presented  in  Sections  4.3  -  4.5  have  been  formatted  for  easier  read¬ 
ability  in  several  fashions  which  need  to  be  mentioned.  The  first  visual  change  is  that 
the  aircraft  and  its  associated  wing-tip  paths  are  scaled  anywhere  from  1  to  20  times 
the  size  of  an  actual  F-16.  This  was  done  to  enable  the  reader  to  visualize  the  attitude 
of  the  aircraft  along  the  associated  trajectory  and  does  not  affect  the  actual  position 
or  attitude  data  from  which  the  plots  were  created. 

Additionally,  each  of  the  optimal  trajectories  are  presented  through  the  use  of 
three  2-D  plots  and  one  3-D  plot.  Through  this  combination  it  is  possible  for  the 
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Aileron  (deg)  Throttle  (pet)  Aileron  (deg)  Throttle  (pet)  Aileron  (deg)  Throttle  (pet) 


(a)  Controls  at  30k  Iterations 


(b)  Trajectory  at  30k  Iterations 


(c)  Controls  at  100k  Iterations 


(e)  Controls  at  146k  Iterations 


(d)  Trajectory  at  100k  Iterations 


(f)  Trajectory  at  146k  Iterations 


Figure  4.3:  Stable  optimization  trajectory  development  over  time. 


56 


reader  to  glean  a  good  understanding  of  each  of  the  aircraft  state  variables  with 
the  exception  of  the  engine  lag  state,  pow.  The  full  state  and  control  time  histories 
associated  with  each  of  the  trajectories  are  located  in  Appendix  A  for  more  detailed 
investigations  of  what  is  occurring  along  the  trajectory. 

Finally,  at  this  point  the  three  aircraft  models,  as  revisited  again  in  Table  4.1, 
performed  each  maneuver  in  a  very  similar  fashion.  Though  there  were  slight  dif¬ 
ferences  in  each  model’s  respective  trajectory,  the  differences  between  them  were  not 
enough  to  warrant  separate  presentations  of  the  results.  To  that  end,  unless  otherwise 
noted,  each  of  the  trajectories  presented  in  this  paper  depict  the  optimal  trajectory 
for  the  baseline  F-16  aircraft  otherwise  known  as  Model  //  1. 


Table  4.1:  Summary  of  Aircraft  Models  Deviation  from  Baseline  Model 


Parameter 

Model  #  1 

Model  #  2 

Model  #  3 

Weight  (W) 

20,500 

20,500 

20,500 

lbs 

Max  Thrust  (T) 

20,000 

25,000 

20,000 

lbs 

Wing  Area  (S') 

300 

300 

375 

ft2 

Wing  Loading  (W / S') 

68.33 

68.33 

54.67 

lbs/  ft2 

Thrust-to- Weight  Ratio  (T/W) 

0.9756 

1.2195 

0.9756 

- 

4-3  Demonstration  Maneuvers 

4-3.1  Northing  Position  Change.  The  results  for  the  Northing  Position 
Change  maneuver,  as  described  in  Table  3.10  and  seen  in  Figure  4.4,  definitely  repre¬ 
sent  the  solution  to  the  classic  brachistochrone  problem.  For  this  maneuver,  the  air¬ 
craft  first  pitches  nose  down  and  trades  altitude  for  increased  velocity  before  pulling 
out  level  to  finish  with  a  dash  to  the  target  state.  Aside  from  the  initial  portion  of  the 
maneuver,  the  resulting  trajectory  follows  a  path  very  similar  to  the  cycloid  solution 
of  the  classic  brachistochrone  maneuver. 

Even  though  the  trajectories  of  the  three  aircraft  were  very  similar  in  basic 
shape,  the  resulting  minimum  times  did  differ  slightly.  The  time  required  to  accom¬ 
plish  the  maneuver  for  each  aircraft  is  found  in  Table  4.2.  Along  with  the  basic  time 
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Figure  4.4:  Northing  position  change  results.  (Model  #  1) 


measurement,  the  final  time  is  also  normalized  against  the  final  time  of  the  Model 
#  1  aircraft  for  comparison  purposes.  These  results  show,  as  expected,  that  the  in¬ 
creased  thrust  of  Model  #  2  allows  it  to  cover  the  distance  faster  and  the  increased 
drag  caused  by  the  wing  area  increase  of  Model  #  3  make  it  the  slowest  of  the  three 
models  for  this  maneuver. 


Table  4.2:  Northing  position  change  maneuver  times. 

Aircraft  Final  Time  Normalized  Time 


Model  #  1 
Model  #  2 
Model  #  3 


15.276  s 
14.747  s 
15.534  s 


1.000 

0.965 

1.017 


4-3.2  Unconstrained  3-D  Position  Change.  The  trajectory  resulting  from 
the  unconstrained  3-D  position  change  maneuver,  as  described  in  Table  3.11  and  seen 
in  Figure  4.5,  resembles  a  classic  Low  Yo-Yo  maneuver.  The  Low  Yo-Yo,  as  perfected 
by  Chinese  fighter  pilot  Yo-Yo  Noritake,  is  an  out-of-plane  lead-pursuit  maneuver 
where  a  pilot  pulls  the  nose  of  the  aircraft  down  and  inside  of  the  turn  in  an  attempt 
to  both  tighten  the  turn  and  increase  speed  at  the  same  time  [24],  In  this  specific 
case,  the  aircraft  rolls  past  100°  of  bank  while  also  dropping  the  nose  as  much  as  15° 
below  the  horizon.  Subsequently,  the  aircraft  slowly  returns  to  wings  level  as  it  pulls 
out  of  the  dive  on  a  course  towards  the  final  target  position. 


Easting  Position  (ft) 

(a)  East-North  Plane  (b)  East-North-Up  View 


(c)  East-Up  Plane  (d)  North-Up  Plane 

Figure  4.5:  Unconstrained  3-D  position  change  Results.  (Model  #  1) 
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Once  again,  though  the  differences  in  final  times  between  the  three  aircraft  are 
very  small,  the  higher  thrust  of  Model  #  2  allows  it  to  complete  the  maneuver  slightly 
faster  than  the  other  two,  as  depicted  in  Table  4.3. 

Table  4.3:  Unconstrained  3-D  position  change  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  21.994  s  1.000 

Model  #  2  20.721  s  0.942 

Model  #  3  21.619  s  0.983 

4-3.3  Constrained  3-D  Position  Change.  Building  on  the  previous  maneu¬ 
ver,  each  of  the  aircraft  accomplish  the  3-D  position  change  maneuver  by  similarly 
rolling  and  dropping  the  nose  of  the  aircraft  to  perform  a  tighter  turn  and  trading 
altitude  for  speed.  The  main  difference  in  this  maneuver,  as  described  in  Table  3.12 
and  seen  in  Figure  4.6,  is  the  rapid  bank  angle  and  heading  change  at  the  end  of  the 
maneuver  to  allow  the  aircraft  to  return  to  its  initial  state  minus  the  position  values. 

As  with  the  previous  two  maneuvers,  the  increased  thrust  available  to  Model 
#  2  enables  it  to  beat  the  other  two  aircraft  to  the  final  position  in  this  maneuver, 
as  found  in  Table  4.4.  Unlike  the  previous  maneuvers,  the  detrimental  effect  of  the 
increased  drag  of  Model  ^  3  is  offset  by  the  increased  maneuverability  which  allow  it 
to  best  the  baseline  model  as  well. 

Table  4.4:  Constrained  3-D  position  change  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  23.891  s  1.000 

Model  #  2  22.962  s  0.961 

Model  #  3  23.396  s  0.979 

4-4  Agility  Maneuvers 

4-4-1  Bank  Angle  Capture.  Kicking  off  the  agility  maneuver  category,  the 
bank  angle  capture  maneuver  is  not  simply  a  90°  roll  as  may  be  deduced  from  the 
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Easting  Position  (ft) 


(a)  East-North  Plane  (b)  East-North-Up  View 


(c)  East-Up  Plane  (d)  North-Up  Plane 

Figure  4.6:  Constrained  3-D  position  change  results.  (Model  #  1) 

trajectories  in  Figure  4.7,  but  a  much  more  complicated  maneuver  due  to  the  other 
target  state  requirements  depicted  in  Table  3.13. 

As  expected,  the  increased  wing  area  and  the  associated  effective  control  surface 
area  increase  allows  Model  ^  3  to  accomplish  this  maneuver  much  faster  than  the 
other  two  aircraft,  as  seen  in  Table  4.5.  The  Model  #  2  aircraft  pulls  out  the  second 
best  time  by  achieving  a  slightly  higher  velocity,  even  in  the  short  maneuver  times, 
which  creates  a  higher  dynamic  pressure  and  therefore  better  control  performance. 
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Figure  4.7:  Constrained  bank  angle  capture  results.  (Model  #  1) 


Table  4.5:  Constrained  bank  angle  capture  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  0.835  s  1.000 

Model  #  2  0.835  s  0.999 

Model  #  3  0.754  s  0.903 


4.4.2  Unconstrained  Heading  Capture.  The  results  for  this  maneuver  simi¬ 
larly  favored  Model  7^  3  with  its  increased  control  authority,  which  is  caused  by  the 
wing  area  being  a  scaling  factor  for  the  control  effect  coefficients.  The  resulting  tra¬ 
jectory  can  be  described  as  a  rapid  roll  and  pull  with  the  aircraft  achieving  the  final 
heading  angle  while  maintaining  a  high  angle  of  attack,  and  still  traveling  along  a 
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Altitude  (ft) 


path  which  deviates  only  slightly  from  its  original  heading.  Figure  4.8  depicts  the 
time-optimal  trajectory  and  Table  4.6  lists  the  corresponding  final  times  associated 
with  this  maneuver  as  numerically  defined  in  Table  3.14. 
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Figure  4.8:  Unconstrained  Heading  Capture  Results.  (Model  #  1) 

Table  4.6:  Unconstrained  heading  capture  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  1.657  s  1.000 

Model  #  2  1.595  s  0.963 

Model  #  3  1.453  s  0.877 

In  an  attempt  to  demonstrate  other  potential  uses  for  this  type  of  trajectory 
optimization  system,  two  addition  variations  of  the  Unconstrained  Heading  Capture 
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were  investigated.  The  first  variation  was  created  by  varying  the  target  heading  angle 
where  —180°  <  ij>  <  180°.  The  results  of  this  investigation  are  found  in  Figure  4.9 
(a)  where  the  radial  distance  from  the  center  depicts  the  time  required  to  complete 
an  Unconstrained  Heading  Capture  of  that  specific  heading  angle.  As  expected,  these 
results  show  that  larger  changes  in  the  heading  angle  result  in  increased  final  times. 


One  interesting  side  note  from  these  results  is  the  asymmetry  in  the  result  data. 
Upon  further  investigation,  it  was  found  that  this  asymmetry  is  the  direct  result  of 
asymmetries  in  the  aerodynamic  model  data.  Two  examples  of  this  phenomenon  are 
found  in  Figure  4.10.  Figure  4.10  (a)  depicts  the  yawing  moment  coefficient,  Cn  as 
a  function  of  the  angle-of-attack,  a,  and  the  sideslip  angle,  /3,  and  Figure  4.10  (b) 
depicts  the  rolling  moment  cause  by  an  aileron  input,  Si/5a,  also  as  a  function  of 
a  and  (3.  Note  that  the  data  is  definitely  not  symmetric  about  (3  =  0  as  would  be 
expected  for  a  symmetric  aircraft. 


(a)  Target  Heading  Variation 
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Figure  4.9:  Variation  of  the  Unconstrained  Heading  Capture  Maneuver. 


The  second  variation  on  the  45°  heading  change  was  to  vary  the  maneuver’s 
initial  conditions.  Figure  4.9  depicts  the  results  of  varying  the  initial  velocity  from 
300  ft/s  to  800  ft/s  and  the  initial  altitude  from  Sea  Level  to  35,000  ft.  With  the 
vertical  axis  representing  the  resulting  final  time,  the  results  show  that,  as  expected, 
this  maneuver  can  be  accomplished  faster  at  higher  initial  velocities  and  lower  initial 
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altitudes.  The  slowest  trajectory  was  found  at  the  point  where  the  aircraft  started  at 
300  ft/s  and  35,000  ft  initial  altitude,  at  which  the  aircraft  is  actually  very  close  to 
the  stall  conditions. 
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(a)  Cn  Asymmetry 


(b)  Si/8a  Asymmetry 


Figure  4.10:  Aerodynamic  Model  Asymmetries. 


4-4-3  Constrained  Heading  Capture.  The  constrained  heading  capture  ma¬ 
neuver  required  the  aircraft  to  both  achieve  a  desired  heading  angle  and  return  to 
steady  level  flight.  In  that  capacity,  the  resulting  trajectory  resembles  a  tight  sus¬ 
tained  turn,  as  described  in  Table  3.15  and  seen  in  Figure  4.11.  As  would  be  expected 
for  a  sustained  maneuver  such  as  this,  the  increased  thrust  of  Model  #  2  allows  that 
aircraft  to  maintain  its  energy  level  throughout  the  turn  and  complete  the  maneuver 
before  either  of  the  other  two  aircraft,  as  depicted  in  Table  4.7.  One  final  note  in  this 
table  is  that  there  is  not  a  final  time  listed  for  Model  #  2.  This  is  due  to  the  fact 
that  this  run  did  not  complete  successfully  during  this  data  collection  series.  This 
is  one  of  the  model  and  maneuver  combinations  which  have  necessitated  the  revised 
cost  function. 

4-4-4  Position  Free  Heading  Reversal.  The  general  results  from  the  position 
free  heading  reversal  maneuver  ended  up  matching  the  results  and  predictions  made 
by  Bocvarov  [4]  almost  exactly.  The  basic  trajectory,  defined  in  Table  3.16  and  seen  in 
Figure  4.12,  is  a  classic  Split-S  maneuver  in  which  the  aircraft  rolls  inverted  and  pulls 
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(a)  East-North  Plane  (b)  East-North-Up  View 
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Figure  4.11:  Constrained  heading  capture  results.  (Model  7^  1) 

Table  4.7:  Constrained  heading  capture  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  10.121  s  1.000 

Model  #  2  9.759  s  0.964 

Model  #  3  9.189  s  0.908 

through  in  a  downward  semicircle  until  achieving  the  target  state  [24] .  This  maneuver 
matches  the  solution  found  by  Bocvarov’s  optimization  of  an  F/A-18  aircraft  model. 

Additionally,  in  his  analysis,  Bocvarov  discussed  the  findings  that  the  most 
important  factor  in  performing  a  rapid  reorientation  of  this  nature  was  the  lift  force 
acting  on  the  aircraft  [4],  Though  not  directly  addressed  by  Bocvarov,  a  simple 
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extension  of  this  analysis  would  lead  to  the  observation  that  an  aircraft  with  a  lower 
wing  loading  should  be  capable  of  achieving  large  scale  reorientation  maneuvers  faster 
than  an  aircraft  with  a  higher  wing  loading.  This  is  exactly  what  is  seen  in  the  results 
for  this  maneuver  found  in  Table  4.8.  The  lower  wing  loading  of  Model  #  3  allows  it 
to  complete  this  maneuver  almost  13%  faster  than  the  baseline  F-16  model  and  over 
9%  faster  than  Model  #  2  with  its  higher  thrust  to  weight  ratio. 

One  item  of  note  with  this  maneuver  is  its  lateral  deviation  from  the  vertical 
plane  as  seen  in  Figure  4.12(a)  and  (c).  This  phenomenon  is  entirely  based  on  the 
external  limits  imposed  in  the  optimization  scheme’s  box  constraints.  In  this  setup, 
the  aircraft’s  pitch  angle,  9,  is  limited  to  6  —  ±86°  to  avoid  the  singularity  at  9  = 
±90°.  The  lateral  deviation  from  the  vertical  plane,  seen  in  Figure  4.12  as  a  600  ft 
progression  in  the  Easting  Position,  is  caused  by  the  optimization  program’s  inability 
to  push  the  aircraft  to  higher  pitch  angles  without  violating  the  box  constraints. 

Table  4.8:  Position  free  heading  reversal  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  13.623  s  1.000 

Model  #  2  13.080  s  0.960 

Model  #  3  11.854  s  0.870 

4.4.5  Position  Fixed  Heading  Reversal.  Again  corroborating  the  predictions 
made  by  Bocvarov,  the  results  for  the  position  fixed  heading  reversal  are  drastically 
different  from  the  position  free  maneuver.  I11  this  maneuver,  defined  in  Table  3.17 
and  found  in  Figure  4.13,  the  aircraft  performs  a  maneuver  which  can  be  likened 
to  a  classic  Hi  Yo-Yo  maneuver  [24],  The  maneuver  commences  with  the  aircraft 
performing  a  small  negative  heading  change  before  pulling  vertical  to  trade  kinetic 
energy  for  potential  energy  as  it  attempts  to  stop  its  forward  velocity.  Inverted  and 
approaching  stall  near  the  top  of  the  arc,  the  aircraft  has  successfully  reversed  its 
velocity  vector  and  continues  to  regain  speed  and  energy  as  it  approaches  the  bottom 
of  the  loop  and  the  terminal  state. 
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Figure  4.12:  Position  Free  Heading  Reversal  Results.  (Model  #  1) 


One  interesting  observation,  also  predicted  by  Bocvarov,  is  that  the  optimal 
trajectory  calls  for  maximum  throttle  for  the  majority  of  the  maneuver.  Though  this 
leads  to  the  fact  that  Model  #  2  can  perform  the  maneuver  faster  than  the  baseline 
F-16,  the  superior  lift  characteristics  of  Model  #  3  again  allow  it  to  perform  this 
maneuver  nearly  6%  faster  than  Model  #  2,  as  seen  in  Table  4.9. 


Table  4.9:  Position  fixed  heading  reversal  times. 

Aircraft  Final  Time  Normalized  Time 

37.352  s 
35.049  s 
32.580  s 


Model  #  1 
Model  #  2 
Model  #  3 
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Figure  4.13:  Position  Fixed  Heading  Reversal  Results.  (Model  1) 


4-4-6  Position  Free  Heading  Reversal  with  Altitude  Floor.  The  addition  of 
a  altitude  floor  restriction  changes  the  position  free  heading  reversal  from  an  out-of¬ 
plane  Split-S  maneuver  to  an  in-plane  sustained  turn,  as  seen  in  Figure  4.14.  With 
the  floor  restrictions,  the  advantage  in  this  maneuver  shifts  back  to  Model  #  2  with 
its  increased  thrust  and  ability  to  sustain  its  energy  through  tighter  turns,  as  found 
in  Table  4.10. 


4-4-1  Initial  State  Capture.  The  results  for  the  Initial  State  Capture  ma¬ 
neuver,  as  defined  in  Table  3.18,  are  found  in  Figure  4.15.  In  this  maneuver,  the 
aircraft  rolls  inverted  and  pulls  through  a  full  negative  loop  while  trading  altitude  to 
gain  velocity  before  pulling  vertical  again  and  completing  the  loop  as  it  returns  to 
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Figure  4.14:  Position  free  heading  reversal  with  altitude  floor  results.  (Model  #  1) 

Table  4.10:  Position  free  heading  reversal  with  altitude  floor  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  19.599  s  1.000 

Model  #  2  18.547  s  0.946 

Model  #  3  17.007  s  0.868 

its  initial  state.  The  resulting  times  from  the  three  different  aircraft  models,  found 
in  Table  4.11,  depict  the  interesting  twist  that  this  maneuver  throws  into  the  agility 
characterization  problem.  Unlike  the  case  of  the  heading  reversal  maneuvers  where 
the  resulting  trajectories  relied  heavily  on  the  attitude  reorientation  capabilities  of  the 
aircraft,  the  results  from  this  maneuver  actually  favor  Model  #  2  with  its  increased 
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ability  in  vertical  maneuvers.  Though  the  increased  lift  available  to  Model  44  3  allows 
it  to  reach  the  0  =  —180°  point  faster,  the  increased  thrust  of  model  44  2  gives  it  the 
edge  in  regaining  the  lost  altitude. 
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Figure  4.15:  Initial  state  capture  results.  (Model  44  1) 


Table  4.11:  Initial  state  capture  times. 

Aircraft  Final  Time  Normalized  Time 

Model  44  1  30.856  s  1.000 

Model  44  2  27.900  s  0.904 

Model  44  3  28.387  s  0.920 
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4-4-8  Initial  State  Capture  with  Altitude  Floor.  As  in  the  case  of  the  po¬ 
sition  free  heading  reversal,  the  addition  of  a  altitude  floor  restriction  changes  the 
initial  state  capture  maneuver  from  an  out-of-plane  vertical  maneuver  to  an  in-plane 
sustained  turn,  as  seen  in  Figure  4.16.  With  the  floor  restrictions,  the  advantage  in 
this  maneuver  shifts  back  to  Model  #  3,  as  found  in  Table  4.12. 
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Figure  4.16:  Initial  state  capture  with  altitude  floor  results.  (Model  #  1) 


Table  4.12:  Initial  state  capture  with  altitude  floor  times. 

Aircraft  Final  Time  Normalized  Time 

36.999  s 
33.667  s 
32.227  s 


Model  #  1 
Model  #  2 
Model  #  3 
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4-5  Compound  Maneuvers 

4-5.1  4-P°int  Position  Change.  The  results  from  the  4- Point  Position 
Change  maneuver,  as  defined  in  Table  3.19,  are  found  in  Figure  4.17.  As  expected, 
the  resulting  trajectory  commences  with  the  aircraft  deviating  from  a  straight  line 
acceleration  with  a  slight  negative  heading  change.  With  bank  angles  approaching 
90°,  the  aircraft  pulls  through  the  first  waypoint  in  a  tight  turn  changing  its  heading 
rapidly  and  setting  itself  up  to  increase  its  velocity  through  the  final  intermediate 
point  on  its  way  to  the  final  state.  Performing  the  maneuver  in  this  manner  allows 
the  aircraft  to  perform  only  minor  heading  changes  and  therefore  maintain  a  higher 
velocity  on  its  approaches  to  the  final  two  points.  The  final  time  results  for  the  various 
models  can  be  found  in  Table  4.13. 

Table  4.13:  4-Point  position  change  times. 

Aircraft  Final  Time  Normalized  Time 

Model  #  1  30.877  s  1.000 

Model  #  2  26.348  s  0.853 

Model  #  3  26.836  s  0.869 


4-6  Potential  as  Control  Method  and  Agility  Prediction  Tool 

Once  the  final  new  cost  function,  Equation  4.5,  was  implemented,  the  approxi¬ 
mate  solutions  provided  by  the  DIDO  optimization  software  almost  always  matched 
the  “truth’’  results  found  from  propagating  the  state  and  controls  histories  as  shown 
in  Section  4. 1.2. 2.  To  that  end,  the  use  of  trajectory  optimization  as  a  tool  for  devel¬ 
oping  control  histories  to  drive  a  simulation,  though  feasible,  is  not  necessary.  The 
same  results  can  now  be  achieved  by  simply  using  the  state  and  control  trajectories 
from  the  optimization  results  for  analysis  and  visualization  purposes. 

If  the  desire  is  truly  to  drive  a  separate  simulation  system,  then  the  optimization 
software  will  be  required  to  interface  with  that  system  and  run  the  optimization 
routine  on  the  full  scale  system.  Any  simplification  of  the  simulator  model  could 
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(a)  East-North  Plane  (b)  East-North-Up  View 


(c)  East-Up  Plane  (d)  North-Up  Plane 

Figure  4.17:  4- Point  position  change  results.  (Model  7^  1) 

results  in  very  poor  results  when  propagating  the  optimum  control  paths  in  a  different 
simulation.  This  issue  is  illustrated  in  Figure  4.18  where  the  control  histories  found 
for  Model  7^  3  to  perform  the  Fixed- Position  Heading  Reversal  maneuver  are  used  to 
propagate  the  solution  for  both  Model  #  1  (a)  and  Model  7^  3  (b)  side  by  side,  with 
poor  results  for  the  Model  #  1  case. 

Another  key  benefit  of  optimizing  the  trajectory  of  a  full  6-DOF  mathematical 
model  is  that  it  now  enables  the  calculations  of  almost  every  proposed  agility  met¬ 
ric  [17].  Since  the  solution  to  a  trajectory  optimization  problem  provides  time  history 
data  for  each  of  the  state  and  control  parameters,  using  a  6-DOF  model  now  pro- 
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Figure  4.18:  Control  interchangeability  example. 


vides  the  necessary  data  for  calculating  every  metric  except  those  based  on  combat 
effectiveness. 

4-7  Overall  Analysis 

As  expected,  the  various  maneuvers  have  shown  that,  though  each  aircraft  has 
definite  advantages  in  certain  conditions,  none  of  these  aircraft  is  superior  throughout 
all  flight  regimes.  Figure  4.19  shows  a  summary  of  the  normalized  optimum  maneuver 
times  for  each  aircraft  and  each  maneuver.  With  lower  values  being  better,  it  was 
found  that  the  aircraft  with  increased  thrust,  Model  #  2,  performed  better  than  the 
other  aircraft  in  maneuvers  such  as  the  Northing  Position  Change  or  Constrained 
Heading  Capture  where  straight  line  acceleration  or  sustained  maneuvers  were  re¬ 
quired.  On  the  other  hand,  the  aircraft  with  an  increased  wing  area,  Model  #  3, 
performs  best  with  instantaneous  reorientation  maneuvers  such  as  the  bank  angle 
capture  and  large  scale  reorientations  like  the  heading  reversal  maneuvers. 

Since  the  various  aircraft  used  in  this  research  are  completely  fictional,  the 
actual  final  times  depicted  in  Figure  4.19  are  of  little  value.  The  real  implication  of 
the  information  in  Figure  4.19  is  that  the  use  of  a  trajectory  optimization  system 
such  as  DIDO  is  a  valid  and  potentially  very  useful  tool  for  predicting  the  differences 


75 


in  aircraft  agility  characteristics  for  various  aircraft.  Additionally,  this  tool  is  robust 
enough  to  handle  all  but  the  most  abstract  of  the  proposed  agility  metrics. 
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Figure  4.19:  Maneuver  final  times  for  the  various  aircraft. 
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V.  Conclusions  and  Recommendations 


5. 1  Conclusions 

The  objective  of  this  research  was  to  develop  a  trajectory  optimization  system 
which  will  alow  a  user  to  investigate  and  compare  the  agility  characteristics  of  various 
aircraft  by  simulating  a  wide  range  of  maneuvers.  As  demonstrated  in  Chapter  4,  a 
system  was  developed  which  allows  a  user  to  simulate  a  wide  variety  of  maneuvers 
with  enough  accuracy  to  render  further  investigations  unnecessary.  In  that  light,  this 
research  was  fundamentally  successful. 

Initial  results  gleaned  from  using  the  baseline  minimum  time  cost  function, 
though  promising,  illustrated  an  instability  in  the  interface  between  the  optimization 
scheme  and  the  mathematical  model.  This  instability  was  characterized  by  increasing 
magnitude  and  frequency  in  the  control  histories  as  a  result  of  the  system  attempting 
to  shave  insignificantly  small  times  off  of  the  final  time.  Since  the  tolerances  for 
defining  an  optimal  solution  in  DIDO  are  not  accessible  to  the  external  user,  it  was 
necessary  to  construct  a  penalizing  term  for  the  cost  function.  Initial  version  of 
the  penalizing  term  focused  on  penalizing  the  actual  use  of  the  controls.  Though 
these  modifications  provided  solutions  in  almost  every  case,  the  time  required  for 
convergence  increased  drastically  as  a  result  of  the  increased  number  of  methods  of 
affecting  the  systems  final  cost.  Results  from  the  use  of  various  penalizing  terms, 
which  instead  focused  on  smoothing  the  control  histories,  depict  faster  convergence 
times  and  better  results  than  were  previously  achieved.  The  final  cost  function  used 
in  this  research  utilized  a  penalty  term  which  was  based  on  the  statistical  variance  of 
the  difference  between  the  actual  control  history  and  a  smoothed  version  of  the  same 
signal.  Through  the  use  of  this  cost  function,  the  optimization  system  was  stabilized 
and  a  wide  variety  of  maneuvers  were  successfully  simulated. 

In  the  area  of  agility  prediction  in  general,  this  research  has  shown  that  it  is 
now  possible  to  optimize  the  trajectory  of  a  full  6-Degree-of-Freedom  mathematical 
aircraft  model.  It  is  no  longer  necessary  to  make  point-mass  approximations  or  any 
other  major  simplifying  assumptions  which  would  cause  the  results  to  be  useless 
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outside  of  academic  circles.  With  few  exceptions,  it  is  now  possible  to  predict  the  full 
gamut  of  metrics  which  have  been  proposed  to  describe  the  agility  characteristics  of 
an  aircraft.  Additionally,  if  the  mathematical  model  is  a  close  enough  representation 
to  the  aircraft  itself,  the  results  from  a  trajectory  optimization  run  can  now  be  used 
as  actual  flight  control  inputs. 

5.2  Recommendations 

The  results  from  this  research,  and  the  subsequent  discussions  and  conclusions, 
lead  to  several  basic  improvements  and  modifications  which  should  be  made  to  the 
actual  trajectory  optimization  system.  The  inherent  instabilities  in  the  optimization 
scheme  are  the  first,  and  most  obvious,  area  in  need  of  future  work.  Further  work 
is  needed  to  characterize  and  mitigate  these  instabilities  by  determining  if  the  insta¬ 
bilities  are  a  result  of  the  mathematical  aircraft  model,  the  optimization  software,  or 
some  combination  of  the  two.  Once  the  optimization  routine  instabilities  have  been 
addressed,  one  could  take  another  look  at  the  ability  of  the  system  to  handle  unstable 
mathematical  models. 

Initially,  it  was  envisioned  that  this  research  would  result  in  the  development 
of  a  Graphical  User  Interface  through  which  a  user  would  have  easy  access  to  all 
of  the  parameters  required  for  the  numerical  maneuver  definition  and  the  associated 
constraints.  This  would  allow  this  tool  to  be  easily  used  for  classroom  projects  and 
demonstrations.  Due  to  time  constraints,  this  portion  of  the  work  was  only  addressed 
at  a  very  basic  level.  A  tool  of  this  nature  could  be  useful  and  would  be  a  good  area 
for  future  work  which  should  not  entail  much  additional  effort. 

The  trajectory  optimization  system  produced  in  this  research  would  be  an  ideal 
foundation  for  a  wide  variety  of  trajectory  optimization  problems.  Relying  heav¬ 
ily  on  the  basic  DIDO  problem  architecture,  the  basic  interface  resulting  from  this 
research  creates  a  highly  intuitive  means  of  setting  up  an  aircraft  trajectory  opti¬ 
mization  problem.  With  only  slight  modifications,  this  basic  structure  can  be  used 
to  simulate  any  mathematical  model  which  can  be  transformed  into  a  self-contained 
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set  of  MATLAB  script  files.  By  doing  so,  many  of  the  aircraft  path  and  trajectory 
optimization  projects  currently  underway  could  benefit  from  this  intuitive  format. 
For  example,  though  the  problem  definitions  in  this  research  did  not  make  use  of  the 
path  constraint  capabilities,  this  system  could  easily  be  used  for  subjects  such  as  UAV 
path  planning  in  urban  environments  or  path  planning  against  RADAR  threats. 

On  a  slightly  different  note,  by  using  the  path  constraint  capabilities  in  this 
system  to  define  the  observed  flight  path  of  an  actual  aircraft,  the  control  and  state 
history  required  to  get  the  aircraft  to  follow  that  flight  path  could  be  backed  out 
of  the  trajectory  optimization  solution.  This  capability  would  be  beneficial  for  both 
accident  investigation  and  reverse  engineering  applications  where  a  certain  subset  of 
the  required  information  is  not  readily  available. 
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Appendix  A.  Expanded  Simulation  Results 
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Figure  A.l:  Northing  position  change  state  and  control  histories.  (Model  #  1) 
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A. 2  Unconstrained  3-D  Position  Change 
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Figure  A. 2:  Unconstrained  3-D  position  change  state  and  control  histories.  (Model 

#1) 
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A. 3  Constrained  3-D  Position  Change 
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Figure  A. 3:  Constrained  3-D  position  change  state  and  control  histories.  (Model 

#1) 
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A. 4  Bank  Angle  Capture 
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Figure  A. 4:  Constrained  bank  angle  capture  state  and  control  histories.  (Model  # 

1) 
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A. 5  Unconstrained  Heading  Capture 


x  104 


(a)  Controls 


(b)  Trajectory 


(c)  Aerodynamic  States 


(d)  Euler  Angles 


(e)  Position 


(f)  Angular  Rates 


Figure  A. 5:  Unconstrained  heading  capture  state  and  control  histories.  (Model  # 

1) 
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A. 6  Constrained  Heading  Capture 
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Figure  A. 6: 


Constrained  heading  capture  state  and  control 
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A. 7  Position  Free  Heading  Reversal 
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Figure  A. 7: 


Position  free  heading  reversal  state  and  control 


histories.  (Model  #  1) 
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A. 8  Position  Fixed  Heading  Reversal 
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Figure  A. 8:  Position  free  heading  reversal  state  and  control  histories.  (Model  A  1) 
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A. 9  Position  Fixed  Heading  Reversal  with  altitude  floor 
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Figure  A. 9:  Position  free 

histories.  (Model  #  1) 
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A.  10  Initial  State  Capture 
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Figure  A.  10: 


Initial  state  capture  state  and  control  histories. 
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A.  11  Initial  State  Capture  with  altitude  floor 
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Figure  A. 11:  Initial  state  capture  with  altitude  floor  state  and  control  histories. 

(Model  #  1) 
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A.  12  4-Point  Position  Change 
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Figure  A.  12: 


4-Point.  Position  Change  state  and  control  histories. 
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